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Abstract 

This  grant  supported  research  into  quiet-flow  supersonic  wind-tunnels,  between  May 
1990  and  December  1994.  Quiet-flow  nozzles  operate  with  laminar  nozzle-wall 
boundary  layers,  in  order  to  provide  low-disturbance  flow  for  studies  of  laminar- 
turbulent  transition  under  conditions  comparable  to  flight.  Major  accomplishments 
include:  (1)  the  design,  fabrication,  and  performance-evaluation  of  a  new  kind  of 
quiet  tunnel,  a  quiet-flow  Ludwieg  tube;  (2)  the  integration  of  pre-existing  codes 
for  nozzle  design,  2D  boundary-layer  computation,  and  transition-estimation  into  a 
single  user-friendly  package  for  quiet-nozzle  design;  and  (3)  the  design  and  prelimi¬ 
nary  evaluation  of  supersonic  nozzles  with  square  cross-section,  as  an  alternative  to 
conventional  quiet-flow  nozzles.  After  a  brief  summary  of  (1),  a  description  of  (2) 
is  presented.  Published  work  describing  (3)  is  then  summarized.  The  report  con¬ 
cludes  with  a  description  of  recent  results  for  the  Tollmien-Schlichting  and  Gortler 
instability  in  one  of  the  square  nozzles  previously  analyzed. 


1  The  Purdue  Quiet-Flow  Ludwieg  Tube 

A  quiet-flow  Ludwieg  tube  has  been  designed,  constructed,  and  tested  at  Purdue 
over  the  past  5  years.  This  is  the  first  low-cost  quiet-flow  facility  constructed  with 
a  design  that  has  the  potential  to  reach  substantial  quiet  Reynolds  numbers  and 
test-section  sizes.  It  is  also  the  first  quiet-flow  facility  with  good  optical  access 
(except  for  the  defunct  JPL  tunnel  [9]  and  the  small  MSU  facility  [7]).  The  facility 
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Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 
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is  currently  in  operation  with  a  3.8  by  4.3  inch  Mach  4  test  section  that  is  quiet 
to  a  length  Reynolds  number  that  exceeds  400,000.  This  Reynolds  number  already 
allows  testing  of  receptivity  and  roughness  effects.  More  importantly,  the  way  is  now 
open  to  development  of  substantially  higher  Reynolds  number  facilities  with  short 
run-times,  at  low  costs.  Current  plans  call  for  the  installation  in  the  Purdue  Ludwieg 
tube  of  the  obsolete  first-generation  nozzle  blocks  for  the  NASA  Langley  Mach  3.5 
pilot  tunnel.  These  blocks,  which  have  a  6  by  10  inch  exit,  should  allow  achieving  a 
quiet  length  Reynolds  number  of  2.7  million  at  30  psia  stagnation  pressure,  according 
to  reference  [4], 

The  design  and  construction  of  the  facility  is  documented  in  references  [13],  [15], 
and  [16].  The  test  results  were  presented  in  reference  [17],  and  have  been  accepted 
for  journal  publication  [14],  More  recent  results  were  presented  in  reference  [19]. 

2  Integrated  Software  for  Design  of  2D  and  Ax- 
isymmetric  Quiet-Flow  Supersonic  Wind  Tun¬ 
nel  Nozzles 

At  the  commencement  of  this  effort,  it  was  expected  that  a  new  nozzle  would  have 
to  be  designed  and  constructed  for  the  Purdue  Ludwieg  tube.  Although  this  has  not 
yet  been  required  after  all,  considerable  effort  has  been  expended  in  developing  an 
integrated,  easy  to  use  quiet-nozzle  design  system  from  the  several  separate  codes 
that  were  originally  in  use  for  this  purpose  at  Langley.  The  following  section  briefly 
describes  and  documents  this  relatively  user-friendly  quiet-nozzle  design  software. 

2.1  Method-of-Characteristics  Nozzle-Design  Code 

The  state  of  the  art  in  the  design  of  supersonic  wind  tunnel  nozzles  involves  the  use 
of  2D  or  axisymmetric  method-of-characteristics  (MOC)  codes  for  determining  the 
nozzle  shapes  that  result  in  uniform  exit  flow.  A  good  introductory  description  of 
the  basic  problems  is  presented  in  sections  15-5  and  16-4  of  reference  [23].  Although 
fully  three-dimensional  MOC  codes  exist,  they  would  have  to  be  iterated  in  order 
to  produce  uniform  flow  at  the  nozzle  exit,  a  basic  requirement  for  wind  tunnel 
nozzles  [12],  Thus,  consideration  is  here  restricted  to  flows  derivable  from  2D  or 
axisymmetric  MOC  solutions. 

The  Sivells  wind-tunnel  nozzle  design  code  was  selected  for  use  in  1990.  Unlike 
the  custom-modified  code  used  by  Chen  for  the  Langley  nozzle  designs  [21],  it  is 
fairly  well  documented  in  reference  [20],  and  the  source  code  is  available.  Although 
it  is  old  FORTRAN-IY  code,  it  is  possible  to  follow  much  of  the  logic  from  the 
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code  itself.  In  addition,  the  Sivells  code  is  capable  of  generating  nozzles  with  the 
radial-flow  regions  that  are  advantageous  from  a  Gortler  stability  viewpoint  [5]. 

The  manual  for  the  code  devotes  8  pages  to  the  description  of  the  input  deck, 
which  allows  for  many  different  combinations  of  nozzle  geometries.  To  simplify  the 
use  of  the  code,  a  FORTRAN  program  was  written  to  automatically  generate  the  in¬ 
put  deck  from  the  answers  to  a  few  simple  questions  that  are  answered  interactively. 
This  program,  SIVINPUT,  which  is  appended,  sets  most  of  the  parameters  in  the 
code.  The  internal  documentation  contained  in  SIVINPUT  in  the  form  of  comments 
will  have  to  suffice  for  documentation  at  the  present  time. 

Numerous  modifications  have  been  made  to  the  Sivells  code  over  the  past  4 
years,  although  most  of  these  are  fairly  minor.  The  mainframe-style  input /output 
structure  was  modified  to  a  friendlier  form,  in  which  the  input  and  output  hies  are 
automatically  generated  by  adding  3  character  suffixes  to  a  user-supplied  8  charac¬ 
ter  rootname.  If  the  rootname  is  here  called  username,  then  SIVINPUT  generates 
username .  inp  which  is  then  read  by  SIVELLS .  A  subroutine  (WRTOBL;  see  appendix) 
was  added  to  write  selected  output  in  a  form  suitable  for  ready  translation  into  Har¬ 
ris  code  input  (hie  username .  bl).  The  WRTOBL  subroutine  also  writes  data  for  the 
conditions  along  the  centerline  of  the  sidewalls  to  the  hie  username .  cl,  and  in¬ 
tegrates  the  tracks  of  mach  lines  originating  on  the  sidewall,  for  determining  the 
width  of  quiet-how  regions  in  2D  nozzles.  The  WRTOBL  subroutine  also  writes  out 
the  derivatives  of  the  wall  contour,  as  generated  in  the  MOC  algorithm,  in  order 
to  obtain  accurate  values  of  the  local  wall  curvature  for  Gortler  computations.  Fi¬ 
nally,  WRTOBL  was  also  modified  to  call  the  Hopkins-Hill  subroutines  developed  and 
documented  in  reference  [2],  in  order  to  determine  the  shape  of  the  bleed  slot  lip 
upstream  of  the  throat.  Although  this  means  that  a  different  algorithm  is  used 
to  determine  the  transonic  throat  shape  upstream  and  downstream  of  the  throat, 
the  two  algorithms  seem  to  merge  smoothly  and  accurately  into  one  another.  This 
agreement  is  not  surprising,  since  both  the  upstream  and  downstream  transonic  flow 
algorithms  are  based  on  near-sonic  perturbation  theory,  so  they  must  agree  very  well 
near  the  throat  where  the  sound  speed  is  nearly  sonic,  at  least  for  large  radius  of 
curvature  wind-tunnel-type  throats.  Error  trapping  code  was  inserted  to  address 
problems  with  user-friendliness,  when  these  were  encountered.  The  common  block 
references  were  made  consistent  to  reduce  compiler  difficulties.  The  array  sizes  al¬ 
lowed  by  the  code  were  also  increased,  to  improve  accuracy.  Minor  changes  were 
also  made  to  SIVELLS  to  allow  restructuring  it  to  structured  Fortran-77  using  the 
commercial  software  package  FOR-STRUCT  (Cobalt  Blue,  Inc.,  Roswell,  Ga.). 

The  only  significant  change  to  the  Sivells  internal  algorithm  was  made  when  it 
was  determined  that  the  code  would  not  generate  internal  streamlines  downstream 
of  the  radial  flow  region,  for  the  nozzle  shapes  desired  (the  original  code  will  only 
generate  internal  streamlines  for  the  special  case  where  ETAD  =  60.)  This  bug  was 
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fixed  by  duplicating  lines  111  and  112  of  the  AXIAL  subroutine  above  line  43  in 
AXIAL. 

Since  a  nozzle  with  10  internal  streamlines  for  square  nozzle  designs  can  be 
generated  on  a  66MHz  486-class  PC  in  less  than  3  minutes,  the  system  is  very 
practical  for  design  work. 

2.2  Sivells  to  Harris  Interface  Code 

The  program  MAKEBLIN  was  written  to  take  the  output  hie  from  the  Sivells  code, 
username. bl,  and  generate  an  input  hie  for  the  Harris  code,  username .bli.  The 
Reynolds  number  scaling  information  required  is  read  in  from  the  auxiliary  hie 
username .  re,  which  contains  the  throat  radius,  total  pressure  and  temperature,  and 
so  on,  and  must  be  hand-generated.  It  uses  various  defaults  to  generate  a  complete 
input  deck  for  the  Harris  code,  including  the  streamwise  grid-point  locations.  It  is 
also  capable  of  generating  Harris  code  input  for  other  inviscid-how  generators.  The 
code  is  heavily  commented,  as  can  be  seen  from  the  listing  in  the  appendix.  The 
highly  automatic  generation  of  the  complex  input  hies  required  for  the  Harris  code 
make  design  studies  relatively  easy  to  carry  out.  This  code  runs  in  seconds  on  the 
66MHz  PC. 

2.3  Harris  Boundary-Layer  Code 

The  Harris  boundary-layer  code  is  documented  with  a  good  user’s  manual  [8]  and  is 
written  in  structured  Fortran.  It  is  a  good  and  standard  finite  difference  code  with 
which  to  compute  2D  and  axisymmetric  boundary  layers.  For  quiet-nozzle  work,  the 
boundary  layers  are  assumed  to  be  laminar  (since  the  nozzle-wall  boundary  layer  is 
only  of  interest  up  to  the  point  where  it  becomes  transitional).  Thus,  the  turbulence 
model  incorporated  in  the  code  is  not  an  issue. 

Again,  minor  modifications  were  made  to  this  code,  to  ease  input/output.  The 
input/output  hies  are  again  automatically  generated  by  appending  suffixes  to  the 
user-supplied  root  name  (e.g.,  username).  The  input  data  is  read  from  username .  bli, 
and  the  standard  output  is  written  to  the  hie  username  .bio.  In  addition,  the 
surface  conditions  written  using  IPRT  commands  are  written  in  tabular  form  to 
username  .prt,  and  the  prohles  written  using  IPRO  commands  are  written  in  tabular 
form  to  username .  pro.  This  allows  rapid  plotting  of  selected  surface  conditions,  and 
eases  translation  of  the  output  into  a  form  suitable  for  the  transition-estimation  code. 
In  particular,  selected  derivatives  of  prohle  quantities  are  written  to  username  .pro, 
as  generated  in  the  program,  so  that  they  can  be  passed  to  the  transition-estimation 
program  with  greater  accuracy.  Error-trapping  code  was  also  added,  as  needed  to 
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ease  use.  The  modified  code  has  been  used  in  required  undergraduate  course  projects 
with  good  success. 

This  code  also  runs  in  a  few  minutes  on  the  66  MHz  PC,  so  again  the  system  is 
convenient  for  design  purposes. 

2.4  Harris  to  e**MALIK  Translation  Code 

The  program  BLTOSTAB  was  written  to  take  the  output  of  the  Harris  code  and 
translate  it  into  the  binary  file  form  used  by  the  e**MALIK  code.  This  rela¬ 
tively  long  code  generates  the  sophisticated  binary  input  file,  username. bfl,  re¬ 
quired  by  the  stability  code  e**MALIK.  It  does  this  by  reading  the  Harris  code 
output  files,  username. prt  and  username . pro,  along  with  the  wall- curvature  data 
saved  in  username .  bl  from  the  Sivells  code,  and  the  Reynolds  number  scaling  data 
saved  in  username. re.  Besides  the  main  binary  output  file,  BLTOSTAB  generates 
username .  cur  for  checking  the  surface  curvature  computations,  username .  gor  for 
printing  Gortler  number  data,  and  username,  sck  for  printing  information  to  check 
the  computations  for  Sivells  input  data.  BLTOSTAB  can  also  read  in  surface  geome¬ 
try  data  that  is  then  differentiated  numerically  to  determine  radii  of  curvature;  this 
feature  was  implemented  to  allow  running  the  Chen  test  case  for  Gortler  instabil¬ 
ity,  described  below.  In  this  case,  the  file  username. usk  is  generated  for  checking 
computations  performed  using  data  not  obtained  from  the  Sivells  code.  A  listing  of 
BLTOSTAB  is  appended.  For  lack  of  resources,  this  heavily  commented  listing  must 
suffice  for  documentation. 

2.5  Transition-Estimation  Code 

The  Langley  program  e**MALIK  was  used  to  perform  e**N  estimates  of  transition 
location  [11,  10].  Although  the  manual,  reference  [11],  is  old  and  imprecise,  it  is 
apparently  the  best  version  in  existence.  The  manual  has  been  supplemented  by 
examination  of  the  source  code,  and  private  communications  with  Mujeeb  Malik, 
Robert  Spall,  and  Scott  Anders.  Because  of  the  uncertainties  of  using  this  complex 
code  with  limited  documentation,  test  cases  were  run  for  both  the  T-S  and  Gortler 
instabilities,  to  confirm  proper  operation. 

The  e**MALIK  code  was  also  modified  slightly,  again  only  to  change  the  in¬ 
put/output  formatting  somewhat,  and  to  add  error  trapping  code.  The  code  reads 
stability  computation  instructions  from  a  file  piped  into  the  executable  with  the 
LInix  <  command,  and  writes  the  general  output  to  a  file  specified  with  the  LInix 
>  command.  It  writes  summary  data  to  a  file  username .  sum;  the  string  username 
is  passed  to  the  code  through  an  additional  read  statement  at  the  bottom  of  the 
input  file.  The  binary-form  input  data  for  the  boundary-layer  profiles  is  read  from 
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username. bfl.  The  code  was  also  modified  to  automatically  loop  through  several 
different  frequencies  or  spanwise  wavenumbers  in  one  execution,  a  property  the  ver¬ 
sion  furnished  to  this  author  in  1991  did  not  have,  although  the  manual  suggested 
that  it  did.  This  is  the  most  CPU-intensive  code  in  this  quiet-nozzle  design  system. 
The  author  has  most  recently  been  executing  it  on  his  department  server,  a  4-CPU 
Sun  Sparcstation  1000  that  has  256  megabytes  of  RAM.  A  half-dozen  different  fre¬ 
quencies  can  be  studied  over  few  dozen  streamwise  stations  in  typical  runs  that 
can  usually  be  carried  out  overnight.  Thus,  although  this  code  is  CPU-intensive,  it 
remains  practical  to  perform  design  studies  with  it  in  a  modern  workstation  envi¬ 
ronment. 

2.5.1  Test  Case  for  T-S  Instability  Computations 

Use  of  the  e**MALIK  code  for  computation  of  T-S  instabilities  was  benchmarked 
back  in  the  summer  of  1990,  using  test  case  6  from  reference  [10].  Test  case  6  is 
for  the  flow  on  a  flat  plate  at  Mach  4.5,  R  =  \J Uex / ue  =  1500,  adiabatic  wall, 
and  a  total  temperature  of  1100  Rankine.  Table  IX  on  p.  407  of  reference  [10]  gives 
spatial  case  results  for  a  nondimensional  frequency  u  =  0.23.  There,  the  eigenvalue  a 
is  given  as  (0.2534081,  -0.0024932)  for  the  most  accurate  computation.  The  author’s 
version  of  e**MALIK  was  first  checked  by  running  it  with  the  internal  self-similar 
boundary-layer  solver,  which  produced  (0.253397,  -0.00250489).  This  agreement  was 
considered  to  be  very  good.  However,  when  the  self-similar  boundary-layer  profiles 
were  generated  by  the  author  of  this  paper,  translated  into  e**MALIK  input  form 
by  BLT0STAB,  and  the  stability  was  recomputed,  the  results  were  not  as  good.  In  this 
case,  at  a  slightly  different  R  =  1495,  a  was  found  to  be  (0.25195,-0.00229).  This  10 
percent  difference  was  doubtless  caused  by  imprecise  generation  of  the  boundary- 
layer  profiles;  the  overall  agreement  indicates  that  the  BLT0STAB  code  generates  the 
proper  e**MALIK  input.  At  the  time  of  this  test,  the  BLT0STAB  code  differentiated 
the  Harris  velocity  profiles  numerically,  instead  of  using  the  internally-generated 
Harris  derivatives,  so  agreement  would  no  doubt  be  better  at  the  present  time. 

2.5.2  Test  Case  for  Gortler  Instability  Computations 

Benchmark  data  for  the  Gortler  test-case  was  generously  supplied  by  Frank  Chen, 
along  with  data  for  the  nozzle  coordinates  and  pressure  distribution.  The  data  is  for 
the  Mach  6  NTC  nozzle  currently  in  use  at  Langley  [6],  at  a  total  pressure  of  100  psia 
and  a  total  temperature  of  360F,  under  adiabatic  wall  conditions.  The  summary 
output  file  printed  by  the  e**MALIK  code  was  supplied,  although  detailed  data 
on  the  boundary-layer  computations  performed  by  Chen  were  no  longer  available. 
This  test-case  tests  both  the  boundary-layer  computations  (carried  out  in  both  cases 
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with  the  Harris  code)  and  the  use  of  the  e**MALIK  code.  Since  the  Chen  nozzle 
design  uses  a  modified  Nelms  code,  and  the  system  described  here  uses  the  Sivells 
code,  the  MOC'  nozzle-design  code  itself  is  not  tested.  Since  the  Chen  data  for  the 
surface  coordinates  are  given  to  four  places,  without  derivatives,  auxiliary  code  was 
added  to  the  BLTOSTAB  routine  to  allow  reading  in  this  non-Sivells  coordinate  data. 
The  Chen  coordinate  data  was  differentiated  using  the  same  LaRC  routine  used  by 
(  hen,  CSDS,  which  was  obtained  from  LaRC  computing  center  personnel. 

f  igure  1  shows  the  radius  of  curvature,  R,  as  a  function  of  the  axial  location 
which  is  zero  at  the  throat.  Note  that  c  is  the  coordinate  along  the  centerline 
of  the  tunnel,  not  the  arclength  s  along  the  nozzle  wall.  Smooth  distributions  are 
obtained  in  the  concave  region  downstream  of  about  2  =  8  inches,  but  the  distribu¬ 
tions  become  noisy  near  the  downstream  end  of  the  nozzle.  Additional  smoothing 
in  the  differentiation  might  have  reduced  this  effect.  The  precise  curvature  distri¬ 
butions  computed  by  Chen  were  unavailable.  The  Harris  code  was  run  using  870 


Figure  1:  Radius  of  Curvature  for  Mach  6  Test  Case 
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streamwise  stations  and  51  points  in  the  wall-normal  direction.  Figure  2  shows  the 
edge  Mach  numbers,  .V/,.  interpolated  by  the  Harris  code  to  the  C'hen-supplied  data. 
The  limited  precision  of  the  data  supplied  causes  the  Harris  code  to  generate  some 


waviness  while  interpolating.  Figure  3  shows  the  boundary- layer  momentum  thick¬ 
ness,  Reg,  comparison.  The  waviness  in  the  Me  interpolation  clearly  causes  some 
waviness  in  the  momentum  thickness  Reynolds  number,  Reg,  and  the  current  results 
are  somewhat  above  those  used  by  Chen.  The  cause  of  the  small  difference  is  diffi¬ 
cult  to  determine,  since  the  details  of  the  Chen  computation  are  no  longer  available. 
Both  the  wall  radius  of  curvature  and  the  momentum  thickness  are  reflected  in  the 
Gortler  number  computations,  presented  in  Figure  4.  Here,  the  Gortler  number  G$ 
is  based  on  the  momentum  thickness  6.  Overall,  the  agreement  is  good,  although 
the  interpolations  and  differentiations  required  in  the  re-analysis  of  the  Chen  data 
clearly  cause  additional  scatter.  No  smoothing  has  been  applied  to  these  computa¬ 
tions.  Finally,  Figure  5  shows  the  results  of  the  N-factor  computations  carried  out 
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Figure  3:  Comparison  of  Ree  for  Mach  6  Test  Case 

using  the  e  MALIK  code,  for  the  case  with  125  waves  around  the  circumference  of 
the  nozzle,  which  is  nearly  the  most  amplified  case.  Here,  s  is  the  arclength  along 
the  nozzle  wall,  with  s  =  0  being  located  at  the  throat.  Clearly,  the  operation  of  the 
Harris,  BLTOSTAB,  and  e**MALIK  codes  is  fundamentally  correct,  for  the  results 
agree  well.  Although  the  close  agreement  is  clearly  in  some  part  fortuitous,  given 
the  minor  disagreements  in  the  previous  plots,  the  test-case  shows  that  this  part  of 
the  integrated  and  fairly  automatic  software  produces  the  proper  results. 

2.6  Summary  of  Software  Status 

This  software  is  now  working  and  fairly  well  tested.  With  this  semi-automatic  soft¬ 
ware,  it  is  possible  to  generate  nozzle  designs  in  seconds  for  2D  and  axisymmetric 
nozzles,  although  square  nozzles  take  longer.  The  boundary  layers  on  2D  and  ax¬ 
isymmetric  nozzles  can  be  computed  in  minutes  on  a  PC,  and  then  a  few  overnight 
computations  on  a  modern  workstation  can  be  performed  to  estimate  the  location 
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z,  inches 

Figure  4:  Comparison  of  G9  for  Mach  6  Test  Case 


of  transition.  Thus,  one  complete  iteration  of  a  quiet  nozzle  design  can  be  carried 
out  in  a  workstation  environment  in  a  few  days,  using  only  an  hour  or  two  per  day 
of  engineering  labor. 


3  Design  and  Preliminary  Evaluation  of  Super¬ 
sonic  Wind  Tunnel  Nozzles  with  Square  Cross 
Sections,  for  Use  in  Quiet-Flow  Facilities 

3.1  Introduction 

This  effort  began  in  September  1992,  and  continued  intermittently  until  the  termi¬ 
nation  of  the  grant.  Sivell’s  code  was  again  used  for  the  method-of-characteristics 
(MOC)  design  of  the  nozzles  [20].  Although  reference  [20]  suggests  that  the  code 
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Figure  5:  Comparison  of  Gortler  N-factors  for  Mach  6  Test  Case 

will  generate  internal  flow  streamlines  for  all  the  types  of  nozzles  which  the  code  can 
generate,  minor  modifications  proved  to  be  necessary  in  order  to  generate  internal 
streamlines  for  nozzles  that  included  radial-flow  regions.  A  post-processing  code  was 
then  written  to  trace  the  streamlines  upstream  from  a  square  cross-section  drawn  in 
the  nozzle-exit  plane,  in  order  to  generate  near-square  cross-sections  from  the  exact 
MOC  results,  between  the  exit  and  the  throat.  This  post-processing  code  was  com¬ 
pleted  in  fall  1992  and  documented  in  reference  [18].  The  Hopkins-Hill  technique  for 
designing  the  transonic  region  of  the  nozzles  (up  to  the  bleed  slot  lip)  still  needed  to 
be  implemented  in  a  reasonably  convenient  way.  This  work  was  completed  during 
early  summer  1993  and  documented  in  reference  [2].  Graduate  student  Timothy 
Alcenius  carried  out  the  Hopkins-Hill  work,  and  then  performed  3D  Navier-Stokes 
computations  of  the  mean  flow  in  the  square  nozzles.  The  position  of  transition  due 
to  the  crossflow  instability  was  estimated  using  these  computations  and  the  cross- 
flow  Reynolds  number  technique.  Most  of  the  Navier-Stokes  results  are  available  in 
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flow  Reynolds  number  technique.  Most  of  the  Navier-Stokes  results  are  available  in 
reference  [1],  which  appeared  in  June  1994.  Complete  results  are  available  in  refer¬ 
ence  [3],  which  will  be  appended  to  the  final  report  for  NAG- 1-1607.  A  summary  of 
the  results  has  been  submitted  for  publication  in  the  AIAA  Journal  of  Aircraft. 

Although  the  above  analyses  suggested  that  crossflow  instability  would  be  the 
dominant  factor  in  transition  of  the  boundary  layers  on  the  square  nozzle  walls,  com¬ 
putations  of  the  Tollmien-Schlichting  instability  and  the  Gortler  instability  were  also 
carried  out  using  the  e**MALIK  code  [10].  Although  the  actual  nozzle-wall  bound¬ 
ary  layer  is  three-dimensional,  these  computations  were  carried  out  using  the  results 
of  a  2D  boundary-layer  computation  of  the  sidewall  boundary  layer  in  the  symmetry 
plane.  Analyses  were  only  carried  out  for  the  long  Mach  2.4  nozzle  described  in  [1], 
since  this  was  thought  to  be  the  worst  case  for  these  instability  modes,  and  since  the 
short  Mach  2.4  nozzle  appeared  less  practical  from  the  standpoint  of  the  crossflow 
instability.  A  detailed  description  of  these  computations  forms  the  remainder  of  this 
section. 

3.2  Computations  of  the  Centerplane  Boundary  Layer  in 
the  Long  Mach  2.4  Square  Nozzle 

The  coordinates  of  the  long  Mach  2.4  nozzle  analyzed  by  Alcenius  are  given  on  p.  129 
of  the  appendix  of  reference  [3].  These  coordinates  were  obtained  by  Alcenius  using 
SIVINPUT  and  the  Sivclls  design  program  [20].  The  input  parameters  are  specified 
on  p.  23  of  reference  [3]  (Nozzle  2  of  Table  1).  It  should  be  noted  that  the  nozzle 
design  codes  produce  a  nozzle  with  a  bleed  slot  lip  that  protrudes  far  upstream  of 
the  throat.  Alcenius  then  cut  off  the  upstream  extent  of  the  bleed  slot,  according  to 
specifications  provided  by  Ivan  Beckwith  from  NASA  Langley,  in  order  to  produce 
a  bleed  slot  originating  0.175  meters  (0.574  ft.)  upstream  of  the  throat  (p.  129  of 
reference  [3]).  For  the  computations  presented  here,  the  bleed  slot  was  cut  off  0.612 
ft.  upstream  of  the  throat;  interpolation  to  the  exact  position  used  by  Alcenius  was 
not  performed.  Figure  6  shows  that  the  coordinates  are  identical,  except  for  the 
difference  in  leading  and  trailing  extent  (Alcenius’s  last  point  is  at  4.197  ft.,  the 
last  point  used  here  is  4.165  ft.).  Here,  z  is  again  the  coordinate  along  the  nozzle 
centerline,  beginning  at  the  throat,  and  the  y-axis  is  normal  to  the  z-axis.  Figure  7 
shows  the  displacement  thickness  (5*)  of  the  2D  boundary  layer  calculated  using  the 
centerplane  pressure  distribution  and  the  Harris  code  [8] ;  the  close  agreement  shows 
that  the  methods  used  were  the  same.  Although  the  small  discrepancy  about  3  ft. 
downstream  of  the  throat  is  somewhat  troubling,  this  slight  variation  in  a  strong 
favorable  pressure  gradient  should  not  have  much  effect  on  the  stability.  Although 
the  point  is  not  discussed  in  reference  [3] ,  it  should  be  noted  that  both  computations 
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Figure  6:  Coordinates  of  Long  Mach  2.4  Square  Nozzle 

computations.  Also,  both  computations  presented  here  were  carried  out  assuming 
an  isothermal  wall  temperature  of  520  Rankine.  Since  this  is  nearly  equal  to  the 
stagnation  temperature,  540  R,  and  the  Mach  number  is  not  high,  the  conditions 
are  nearly  adiabatic.  The  stagnation  pressure  is  taken  as  100  psia,  as  in  Alcenius’s 
work. 

3.3  Estimates  of  Transition  in  the  Long  Mach  2.4  Square 
Nozzle  due  to  the  Gortler  Instability 

This  instability  was  analyzed  with  the  e**MALIK  code,  using  the  nozzle-design 
software  described  above.  The  boundary- layer  data  were  read  in  through  the  exter¬ 
nal  binary  file  generated  from  the  Harris  code  output  by  the  translation  program 
BLTOSTAB .  The  body  was  assumed  2D,  71  grid  points  were  used  in  global  computa¬ 
tions  and  141  in  local  computations,  and  computations  were  begun  at  the  location 
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Figure  7:  Comparison  of  Centerplane  6*  Values  for  Long  Mach  2.4  Square  Nozzle 


where  the  curvature  becomes  concave.  Since  ibeta  was  set  to  0,  the  input  parameter 
betx  is  actually  the  dimensional  wavenumber,  k,  per  ft.  (note  that  the  e**MALIK 
manual  [11]  is  erroneous  in  this  regard;  the  correction  was  provided  verbally  by 
M.  Malik  on  16  July  1991).  This  wavenumber  parameter  was  varied  through  500, 
1000,  1200,  1400,  1600,  1800,  2000,  2200,  2400,  2600,  2800,  3000,  and  4000  per 
ft.  Since  the  momentum  thickness  8  of  the  boundary  layer  at  the  first  concave 
station  was  about  2  x  10-4  ft.,  and  the  wavelength  A  =  2i r/k,  this  corresponds  to 
spanwise  wavelengths  X/8  ranging  from  60  to  8.  The  corresponding  values  of  BETA 
output  by  the  code  range  from  0.19  to  1.9;  these  spanwise  wavenumber  values  are 
non-dimensionalized  by  the  distance  from  the  leading  edge  and  the  edge  values  of 
the  fluid  properties.  The  spanwise  wavelength  of  the  Gortler  waves  is  maintained 
constant  as  the  amplification  is  integrated  downstream.  Figure  8  summarizes  the 
integrated  amplification  between  the  streamwise  distances  of  2.80  and  4.17  ft.  The 
peak  amplification  is  about  5.7,  which  occurs  at  a  BETA  of  0.9.  Based  on  the  Langley 
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Figure  8:  Gortler  N-factors  at  Exit  of  Long  Mach  2.4  Square  Nozzle 


N-factor  criterion  of  7.5  for  Gortler  transition  [22],  this  instability  is  not  expected 
to  cause  transition.  Of  course,  this  estimate  neglects  interactions  with  roughness  or 
other  instabilities. 


3.4  Estimates  of  Transition  due  to  the  T-S  Instability 

The  T-S  computations  were  performed  with  the  e**MALIK  code,  and  with  the  same 
boundary-layer  profiles  used  for  the  Gortler  computations.  Many  computations  had 
to  be  performed  in  order  to  find  the  unstable  range  of  frequencies.  The  results  are 
shown  in  Figure  9.  Frequencies  in  the  neighborhood  of  6  kHz  are  most  unstable. 
These  have  a  wave-angle  PS  I  of  about  67  degrees,  which  appears  to  be  reasonable. 
The  nondimensional  frequency  OMEGA  is  about  0.009  for  these  waves.  Since  the 
integrated  N-factor  is  less  than  2,  it  appears  unlikely  that  T-S  instability  would 
contribute  substantially  to  transition  in  this  nozzle. 
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arclength  s,  ft. 

Figure  9.  T-S  N-Factors  in  Long  Mach  2.4  Square  Nozzle 


3.5  Summary  of  Square-Nozzle  Design  Effort 


The  design,  construction,  and  shakedown  of  high  speed  quiet-flow  nozzles  is  a  dif¬ 
ficult  task  in  which  many  poorly  understood  tradeoffs  have  to  be  made.  All  of  the 
difficult  high-speed  transition  phenomena  that  the  nozzles  are  intended  to  study 
must  be  estimated  in  order  to  best  design  the  nozzles.  Fabrication  costs  are  difficult 
to  estimate  in  advance. 


The  square  nozzle  concept  put  forward  by  Ivan  Beckwith  of  NASA  Langley  was 
suggested  as  a  possible  method  of  obtaining  the  advantages  of  both  axisymmetric 
and  rectangular  (2D)  nozzles.  The  walls  are  flat  downstream  of  the  expansion, 
facilitating  the  use  of  windows.  The  nozzle  can  be  disassembled  to  facilitate  polishing 
and  maintenance  of  the  critical  throat  region.  High  Mach  number  square  nozzles 
can  be  machined  without  excessive  difficulties  with  throat  tolerances,  and  contour 
flaws  do  not  focus  to  the  centerline.  All  four  walls  see  the  same  accelerating  flow, 
so  it  was  hoped  that  crossflow  would  not  be  a  major  problem. 
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facilitating  the  use  of  windows.  The  nozzle  can  be  disassembled  to  facilitate  polishing 
and  maintenance  of  the  critical  throat  region.  High  Mach  number  square  nozzles 
can  be  machined  without  excessive  difficulties  with  throat  tolerances,  and  contour 
flaws  do  not  focus  to  the  centerline.  All  four  walls  see  the  same  accelerating  flow, 
so  it  was  hoped  that  crossflow  would  not  be  a  major  problem. 

The  studies  carried  out  to  date  do  not  support  the  high  hopes  initially  conceived. 
Estimates  based  on  the  crossflow  Reynolds  number  (Alcenius)  suggested  that  cross- 
flow  would  be  a  major  difficulty  in  these  nozzles,  even  in  the  throat.  These  estimates 
have  apparently  been  confirmed  by  unpublished  e**N  estimates  carried  out  by  the 
High  Technology  Corp.  under  the  SBIR  program.  It  remains  possible  that  the 
throat-region  crossflow-transition  problem  identified  by  these  computations  could 
be  ameliorated  by  moving  the  bleed  slot  downstream.  The  cost  of  machining  such  a 
nozzle  remains  a  topic  of  speculation.  It  may  still  be  that  square  nozzles  will  prove 
to  be  more  cost-effective  than  axisymmetric  or  2D  designs.  The  research  to  date 
seems  somewhat  discouraging  but  far  from  conclusive. 

4  Measurements  of  Crossflow  Instability  on  the 
Sidewalls  of  the  LaRC  Mach  3.5  Low  Distur¬ 
bance  Tunnel 

Part  of  the  proposed  crossflow  investigations  involved  measurements  of  the  flat  side- 
wall  boundary  layers  in  the  Mach  3.5  low  disturbance  tunnel  at  Langley.  These 
were  to  be  carried  out  during  summer  1994  by  Christine  Haven,  as  part  of  her 
M.S.  thesis  at  Purdue.  One  and  a  half  days  of  tunnel  access  were  provided  in  the 
middle  of  the  summer,  after  about  a  month  of  preparation.  Although  the  hot  wires, 
traverse  system,  and  controller  software  worked  acceptably,  major  difficulties  with 
electronic  noise  were  encountered.  Although  the  cause  of  these  difficulties  was  later 
determined,  the  problem  could  not  be  solved  in  time  to  allow  obtaining  low-noise 
data.  The  high  noise  level  and  the  limited  amount  of  data  that  could  be  obtained 
has  precluded  any  attempt  to  draw  definite  conclusions  from  this  work.  A  full  report 
will  be  provided  with  the  final  report  on  Langley  grant  NAG- 1-1607,  which  should 
be  forthcoming  in  a  few  months. 


5  Summary 

This  grant  enabled  the  successful  development  of  a  new  kind  of  low  cost  quiet-flow 
wind  tunnel  at  Purdue.  Since  this  wind  tunnel  is  the  basis  for  current  AFOSR- 
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supported  research,  the  outcome  of  the  grant  has  involved  successful  technology- 
transfer  to  AFOSR.  The  grant  also  supported  the  development  of  relatively  efficient 
and  user-friendly  software  for  quiet  nozzle  design,  although  the  full  potential  of 
this  software  has  yet  to  be  put  to  use.  Finally,  the  grant  supported  the  design  and 
preliminary  evaluation  of  the  new  square  nozzle  concept  for  quiet-flow  wind  tunnels. 
Although  the  preliminary  evaluations  of  this  concept  are  somewhat  discouraging, 
further  research  might  yet  show  that  square  nozzles  will  be  an  advantageous  choice 
in  certain  circumstances. 
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A  Source  Code  for  SIVINPUT 


*  PROGRAM  SIVINPUT.  Fortran-77.  For  generating  input  deck  to 

*  Sivells  nozzle  design  code  (AEDC-TR-78-63) . 

* 

*  emailed  back  ‘From  moen  Wed  Nov  25  12:00:09  1992’  after  use  by  him. 

* 

c  this  program  writes  an  input  file  in  the  correct  format  for 
c  sivells  code;  sps  6-27-90 

c  modified  9-30-92  for  the  square  nozzle  work  with  streamlines  sps 

c  modified  to  produce  starting  streamlines  on  nozzle  wall. 

*  mod  sps  11-30-92  to  make  mp=5 

*  mod  11-30  sps  to  allow  choices  of  IN,  XC,  IX 

*  add  header  info.  1-16-95,  this  is  the  current  version  used  by  Alcenius 

*  for  his  MS  thesis  and  square  nozzle  computations  (see  appendix  of 

*  Alcenius  MS  thesis) . 
c 

character*10  title 
character*20  sivfile 
c 

write  (*,*)  ’enter  a  rootname  to  write  sivells  input  file:’ 
read  (*,10)  sivfile 
title  =  sivfile 
ileng  =  index(sivf ile, 1  ’)  -1 

sivfile (ileng+1 : ileng+4)  =  ’.inp’ 
write  (*,*)  ’opening  file-’ , sivfile, ’-for  output’ 
open(unit=2 , f ile=sivf ile , status=’new’ ) 
c 

*  write  (*,*)  ’enter  title  of  run  (10  characters):  ’ 

*  read  (*,20)  title 

write  (*,*)  ’enter  jd  (-1=2D,  0=axisym) :  ’ 

read  (*,*)  jd 

write  (2,30)  title, jd 

*  write  (*,*)  ’enter  sfoa: ’ 

*  read  (*,*)  sfoa 

sfoa=0.  !use  3rd  or  4th  degree  distribution 
gam  =  1.40 
ar  =  1716.563 
zo=l 

c  following  three  used  in  bl  computations,  not  used  here 
ro=l 
visc=l 
vism=l 
c 

xbl=1000.  ! gives  values  at  evenly  spaced  intervals 

write  (2,40)  gam, ar,zo,ro, vise, vism, sfoa, xbl 
c 

write  (*,*)  ’enter  etad,rc,bmach,cmc:  ’ 
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c 


c 


c 


c 

c 

c 

c 

c 

c 


read  (*,*)  etad,rc,bmach,cmc 
xc=0.  Iso  4th  degree  distribution,  change? 
write  (*,*)  ’enter  xc,  in:  ’ 
read  (*,*)  xc,in 
write  (*,*)  ’xc,in=  ’,xc,in 
fmach=0.  Ithis  sets  distribution,  change? 
sf  =  0 .  ! nozzle  throat  radius  =  1.0 

pp  =  0 .  ! coordinates  given  relative  to  throat 

write  (2,50)  etad.rc ,fmach,bmach, cmc , sf ,pp,xc 
write  (*,*)  ’enter  mt ,nt , ix, in,md,nd,nf ,mp, j c , lr ,nx: ’ 
read  (*,*)  mt ,nt , ix, in,md,nd,nf ,mp, jc,lr,nx 
mt=61  Ipts  on  char.  EG,  max  125 
nt=31  Ipts  on  axis  IE,  max  149-LR 
write  (*,*)  ’enter  ix:  ’ 
read  (*,*)  ix 

ix=0  I  is  3rd  deriv  matched?  change? 

in=10  I  use  Mach  no.  distrib  on  BC,  makes  2nd  deriv  match  rad.  flow 

change? 

iq=0  I  calls  for  complete  contour 

md=61  Ipts  on  char.  AB,  max  about  125,  odd 

nd=15  Ipts  on  axis  BC,  max  about  150, 

changed  from  49  to  15  sps  7-2-90 
write  (*,*)  ’enter  -1  for  smoothing,  1  for  no  smoothing:  ’ 
read  (*,*)  ismooth 

nf =ismooth*81  Ipts  on  characteristic  CD.  Neg  calls  for  smoothing 

mp=5  Ipts  on  GA,  conical  section,  if  Fmach  ne  Bmach 

jc=0  I  if  not  0,  used  to  print  intermediate  characteristics 

lr=31  Ipts  on  throat  char.,  -  prints  out  transonic  soln 

nx=13  I  spacing  of  pts  on  axis  upstream,  this  no.  recc. 

mq=0  Ipts  downstream  of  D 

jb=-l  I neg  for  no  BL  computation 

jx=l  Ipos  calls  for  streamlines 

it=0  I  jack  points,  not  used 

write  (2,60)  mt ,nt , ix, in, iq,md,nd,nf ,mp,mq, jb, jx, jc, 

>  it,lr,nx 

if  (ismooth  .eq.  -1)  then 

noup=10  I  smoothing  parameters,  arbitrary 

nodo=10 

npct=90 

write  (2,70)  noup,npct ,nodo 
end  if 

gives  streamline  distribution  that  corresponds 
to  the  half  wall  for  conversion  to  a  square  nozzle, 
note  that  the  number  of  streamlines  requested  will 
be  reduced  by  one  because  Sivells  automatically 
calculates  the  wall  streamline.  Sivells  output 
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c  will  have  the  actual  number  of  streamlines  requested, 
c  (moen  10-92) 
c 

write  (*,*)  ’How  many  streamlines  along  half wall?1 

read  (*,*)  nstream 

nstream=nstream-l 

dx=l . 0/ (float (nstream) *sqrt (2.0)) 

ycnt=l . 0/sqrt (2 . 0) 

do  100  istream  =  0,nstream-l 

etadstr  =  etad*sqrt ( (istream*dx) **2+ycnt**2)  ! see  btm  page  59 
qm  =  sqrt ( (istream*dx) **2+ycnt**2) 
xj  =  1  Hook  for  more  streamlines 
write  (2,90)  etadstr ,qm,xj 
if  (ismooth  .eq.  -1)  then 

write  (2,70)  noup,npct ,nodo  !must  have  for  each!! 
end  if 
100  continue 
c 

close (2) 
stop 

10  format (a20) 

20  format (alO) 

30  format(lx,al0,2x,i2) 

40  format(8(lx,f9.3)) 

50  format(8(lx,f9.3)) 

60  format(lx,i4, 15(i5)) 

70  format(lx,i4,2(i5)) 

90  format(3(lx,f9.4) ,lx) 
end 

B  Output  Subroutine  Written  for  Sivells  Code 

The  following  is  the  subroutine  added  to  the  Sivcll’s  code  to  generate  output  files 
for  transfer  to  the  boundary-layer  code. 

*  subroutine  wrtobl.for 

*  this  subroutine  modified  from  perfc  by  s  schneider  6-90  to  write  sivells 

*  output  to  a  file  to  be  read  directly  by  a  bl  program. 

*  also  computes  the  upstream  contour  from  halls  assumptions  and  rc 

*  assumes  throat  is  0,0,  and  uses  formulas  from  Hall,  (3)  and 

*  assumes  sf=0.0  at  beginning  so  throat  radius  is  1  unit 

*  modified  7-31-91  to  give  full  reference  to  block  C0NTR  same  as  elsewhere 

*  modified  9-29-92  to  output  multiple  calls  when  writing  streamlines, 

*  also  start  z  write  at  Mach  >  1 

*  dimension  arrays  in  parameter  statement.  Note  arrays  must  match  with 

*  other  routines.  Change  small  from  0.05  to  0.01 
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*  modified  9-30-92  to  find  streamline  shapes  in  transonic  region 

*  mod  11-25-92  to  increase  array  sizes  and  use  parameter  statements  sps 

*  mod  6-4-93  to  call  hophill  to  do  hopkins-hill  upstream  of  throat,  sps 

*  mod  9-1-94  to  write  out  WALTAN  and  sd  with  contour  for  Gortler,  sps 

SUBROUTINE  wrtobl  1 

C  4 

parameter  (mwr=300 ,maxtpt=150 ,mpt=500 ,mwk=400) 

IMPLICIT  REAL*8 (A-H , 0-Z)  5 

character*20  blf ile , elf ile 
logical  lopen 
dimension  pepO(maxtpt) 

*  following  single  for  passing  to  hophill 

real*4  rap, thetap, gamp, xendp.ytp 

real*4  xp(maxtpt) ,yp(maxtpt) ,amachp(maxtpt)  Isingle  for  passing 

*  (upstream  extension  of  contour  into  subsonic-transonic  region) 

save  icall , xend , numpt s , rap , thetap , gamp , xendp 
common  /centerl/clx(mwr) ,clmach(mwr) 


COMMON  /GG/  GAM,GM,G1,G2,G3,G4,G5,G6,G7,G8,G9,GA,RGA,QT  6 

COMMON  /CLINE/  AXIS (5 ,mwk) ,TAXI (5 ,mwk) , WIP ,X1 , FRIP , ZONK, SEO , CSE  7 

COMMON  /COORD/  S(mpt) ,FS(mpt) , WALTAN (mpt) ,SD(mpt) ,WMN(mpt) ,TTR(mpt  8 

1) , DMDX (mpt ) , SPR (mpt ) ,DPX(mpt) ,SECD(mpt) , XBIN , XCIN , GMA , GMB , GMC , GMD  9 

*  following  not  used  here,  messy  common  block,  avoid  if  possible 

*  COMMON  /WORK/  A(5 , 150) , B (5 , 150) , FINAL (5 , 150) , WALL (5 ,mpt) , WAX (mpt ) ,  10 

*  lWAY(mpt) , WAN (mpt)  11 

COMMON  /PROP/  AR, Z0 , R0 , VISC , VISM, SF0A , XBL , C0NV  12 

COMMON  /PARAM/  ETAD , RC , AMACH , BMACH , CMACH , EMACH , GMACH , FRC , SF , WWO , WW  13 

10P ,  QM ,  WE, CBET ,XE, ETA, EPSI  ,BPSI , XO ,  YO , RRC ,  SDO , XB , XC ,  AH , PP ,  SE , TYE , XA  14 

*  COMMON  /TROAT/  FC(6,51)  15 

*  COMMON  /CONTR/  ITLE(3) , IE, LR, IT , JB , JQ , JX , KAT, KBL , KING ,K0 , LV , NOCON ,  16 

*  1IN, MC ,MCP , IP , IQ , ISE, JC ,M,MP ,MQ , N , NP , NF ,NUT , nr , lc ,md ,mf ,mt ,nd,nt  17 

*  this  is  full  common  block  reference  to  CONTR  taken  from  routine  AXIAL 


COMMON  /CONTR/  ITLE(3) , IE, LR, IT , JB , JQ , JX , KAT, KBL , KING , KO , LV , NOCON, AXI  12 
1  IN ,  MC , MCP ,  IP ,  IQ ,  ISE ,  JC , M , MP ,  MQ , N ,  NP ,  NF , NUT ,  NR ,  LC , MD ,  MF , MT ,  ND ,  NT  AXI  13 

* 

*  following  common  block  filled  by  call  to  tcoeff,  gives  coefficients  of 

*  transonic  series  solution  -  sps 

common  /transc/ gr , gs , gt , gv , gk , u42 , u22 , u63 , u43 , u23 , up2 , upO , 

>  v42,v22,v02,v63,v43,v23,v03 

data  tiny/1 . Oe-5/ , eangle/10 . 0/  ! degrees,  must  be  gt .  0 

*  (allow  1  percent  error  in  computation  of  upstream  contour;  this  is 

*  a  small  size  of  mach  number....) 

*  Note  that  a  choice  of  0.05  for  small  gives  amstar2  >  possible  in  upstream 

*  part  of  transonic  solution  near  centerline,  9-29-92 

*  change  small  from  0.01  to  0.002  10-21-92,  was  giving  problems  when 

*  Mike  Moen  was  running  mach  3.5  test  cases  in  axisym 

data  icall/1/ 

C 

if  (icall  .eq.  1)  then 
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*  (first  call  to  wrtobl) 

inquire  (unit=3 , name=blf ile , opened=lopen) 
if  (lopen)  then 

write  (*,*)  ’writing  nozzle  summary  for  bl  read  to  ’, bifile 
else 

write  (*,*)  ’no  file  open  for  writing  bl  info  to!!!???’ 
return 
end  if 

inquire  (unit=4 , name=clf ile , opened=lopen) 
if  (lopen)  then 

write  (*,*)  ’writing  nozzle  summary  for  cl  to  ’,clfile 
else 

write  (*,*)  ’no  file  open  for  writing  cl  info  to!!!???’ 
return 
end  if 
else 

write  (*,*)  ’wrtobl,  call  ’,icall,’  for  streamline  dump’ 
end  if 

* 

*  get  some  things  needed  for  transonic  case  and  for  cl  mach  distribut 

if  (ie  .eq.  0)  then 

*  (planar) 

sigma  =  0.0 
else 

*  (axisymmetric) 

sigma  =  1.0 
end  if 
si  =  rc  +  1 

alambd  =  sqrt (  (l+sigma)/(  (gam+l)*sl  )  ) 

* 

*  write  out  the  cl  mach  number  distribution  info,  for  first  call  only 

*  subsequent  calls  are  to  write  streamline  data 

* 

if  (icall  .eq.  1)  then 

*  first,  compute  the  number  of  points  already  in  the  array: 

*  (this  is  necessary  because  of  overlap  and  half-filling  in  method) 

numpts  =  1 
do  10  i  =  2,mwr 

if  (clx(i)  .gt.  clx(numpts))  then 

*  (a  real  point) 
numpts  =  numpts  +  1 
clx(numpts)  =  clx(i) 
clmach (numpts)  =  clmach(i) 

else 

*  (should  be  a  zero  point) 

if  (clx(i)  .ne.  0.)  then 

write  (2,*)  ’WRTOBL:  dropping  cl  point  which  is  ’, 

i,clx(i) ,clmach(i) 


> 


end  if 
end  if 
10  continue 

write  (4,95)  ITLE,XBIN,XCIN,SF,frip 

write  (4,*)  ’z  computation  begins  at  first  point  where  M>1’ 
write  (4,*)  ’  following  line  is  npts,  then  x.machno. ,z,mu  ’ 
write  (4,*)  numpts  +  2*nt 

*  first,  write  the  points  up  to  the  radial  flow  region: 

z  =  0.0 

do  20  il  =  1, numpts 

if  (clmach(il)  .le.  emach)  then 
if  (clmach(il)  .gt.  1.0)  then 
amu  =  dasin(l . 0/clmach(il) ) 
else 

c  write  (*,*)  ’mach  number  on  centerline  It  1,  skip  when1 

c  write  (*,*)  ’  computing  sidewall  mach  lines’ 

c  write  (*,*)  ’for  clmach,il=’ ,clmach(il) ,il 

amu  =  0.0 

end  if 

if  (il  .gt.  1)  then 

z  =  z  +  tan(amu)*(clx(il)-clx(il-l)) 
end  if 

write  (4,15)  clx(il) ,clmach(il) ,z,amu 
15  format(4(el4.7,2x)) 

illast  =  il 
else 

go  to  21 

*  *(exit  the  loop)* 
end  if 

20  continue 

21  continue 

*  now,  write  the  radial  flow  region 

if  (abs(xbin-(xb*sf+frip))  .gt.  tiny)  then 
write  (*,*)  ’  xbin,xb=  ’,xbin,xb 
write  (*,*)  ’  sf ,frip=  ’,sf,frip 
stop  ’problem  with  xbin’ 
end  if 

xein  =  xe*sf+frip 

c  write  (*,*)  ’xbin,xein=  ’, xbin, xein 

deltam  =  (bmach  -  emach) /(2*nt) 
garni  =  2.0/(gam+l) 
gam2  =  (gam-l)/(gam+l) 
gam3  =  (gam+l)/(2.0*(gam-l)) 
c  write  (*,*)  ’gam, 1 , 2 ,3=’ ,gam,gaml ,gam2 ,gam3 

do  30  i2  =  l,2*nt 

xmach  =  emach  +  i2*deltam 

*  (following  implements  eqn  29  of  sivells  report) 

rhs  =  ((garni  +  gam2*xmach**2) **gam3) /xmach 
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xsorce  =  rhs** (1 . 0/(1 . 0+sigma) ) 
xin  =  xsorce*sf  +  frip 
amu  =  dasin(l . 0/xmach) 
if  (i2  .gt.  1)  then 

z  =  z  +  tan(amu) * (xin-xinold) 
xinold  =  xin 
else 

*  write  (4,*)  ’z,amu,clx(illast) ,xin=’ ,z,amu,clx(illast) ,xin 
z  =  z  +  tan  (amu)  *  (xin-clx(i  Hast)  ) 

xinold  =  xin 
end  if 

write  (4,15)  xin,xmach,z,amu 
30  continue 

c  write  (*,*)  ’sigma=’ , sigma 

c  write  (*,*)  ’ xin, xmach, xsor ce, rhs= ’, xin, xmach, xsorce, rhs 

*  now,  write  out  the  points  beyond  the  radial  flow  region 

do  40  i3  =  il,numpts 

if  (clmach(i3)  .gt.  bmach)  then 
amu  =  dasin(l . 0/clmach(i3) ) 
z  =  z  +  tan(amu)*(clx(i3)-clx(i3-l)) 
write  (4,15)  clx(i3) ,clmach(i3) ,z,amu 
end  if 
40  continue 

close  (unit=4) 

*  (done  with  writing  cl  mach  distribution)  * 

end  if 

* 

*  following  gives  coefficients  of  transonic  soln,  needed  for  pressure 

call  tcoeff(gam) 

* 

if  (pp  .ne.  0)  stop  ’coords  should  be  computed  rel.  to  throat’ 

*  (checks  to  see  if  coords  really  computed  relative  to  throat) 

* 

*  now  compute  the  transonic  extension  upstream: 

*  (note  that  Sivells  array  below  will  begin  at  about  the  throat ! ) 

*  (following  uses  only  that  the  throat  has  a  radius  of  curvature  of 

*  rc  at  the  throat .  Upstream  of  the  throat ,  the  shape  of  the  entrance 

*  is  arbitrary,  as  far  as  the  small  perturbation  transonic  solution  in 

*  the  throat  is  concerned.) 

*  first,  determine  farthest  upstream  can  reasonably  compute: 

*  (let  second  term  in  contour  be  small  compared  to  first) 

*  (this  only  works  for  the  nozzle  contour,  which  is  arbitrary  and  determines 

*  the  interior  streamlines)* 

if  (icall  .eq.  1)  then 

write  (*,*)  ’Using  Hopkins-Hill  nozzle  shape  in  upstream  tr.’ 
write  (*,*)  ’enter  entry  angle,  degrees,  gtO:  ’ 
read  (*,*)  eangle 

write  (*,*)  ’using  nozzle  entry  angle  of  ’, eangle,’  degrees’ 
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*  now  compute  a  series  of  points  down  to  the  throat: 

*  6-28-90  skip  the  point  at  the  throat,  end  up  with  two,  gives  trouble 

*  (compute  same  number  of  points  as  used  in  throat  characteristic) 

numpts  =  abs(lr) 

if  (numpts+1  .gt.  maxtpt)  stop  ’too  many  points  in  wrtobl  array’ 
yth  =  1.0  ! throat  radius  is  1.0,  eangle  is  degrees 
xend  =  -l*tiny 

*  (Ismail  fraction  of  radii,  upstream  near  throat,  not  overlap) 

*  (corrected  sign  9-2-94  sps) 

*  Note  that  hophill  is  single  precision!!  must  convert!! 

rap  =  rc 
thetap=eangle 
gamp  =  gam 
xendp  =  xend 
ytp  =  yth 

call  hophill (rap , thetap , gamp , xendp , ytp , numpts , xp , yp , amachp) 

*  this  routine  returns  numpts  values  of  x,y,amach  along  streamline 

*  upstream  from  eangle  to  xend.  Sets  xbegin  for  later  calls. 

else 

xthroat  =  s(l)  ! these  are  the  first  points  in  Sivells  streamline 
ytp  =  fs(l)  !must  connect  to  transonic  streamline 
if  (xthroat  .ne.  0.0)  then 

write  (*,*)  ’xthroat=  ’.xthroat,’  first  sivells  pt ’ 

*  pause  ’WRTOBL:  x  should  be  0  at  throat’ 

*  first  point  on  streamline  MUST  be  greater  than  0,  for  sonic  line 

*  bows  downstream! 

end  if 

call  hophill (rap , thetap , gamp , xendp , ytp , numpts , xp , yp , amachp) 
end  if 

do  60  itx  =  1, numpts  !get  pepO,  convert  back  to  double 
amach2  =  amachp (itx) **2 

denom  =  (1.  +  amach2* (gam-1 .) /2 . 0) ** (gam/ (gam-1 .) ) 
pepO(itx)  =  1.0/denom 
60  continue 

* 

*  now  write  out  the  bl  data  to  a  special  file  in  easily  read  form 

* 

if  (icall  .eq.  1)  then 

write  (3,95)  ITLE,XBIN,XCIN,SF,frip 
95  FORMAT  (1H1.3A4,’  sivells,  xbin=  ’.F11.6, 

>  ’  xcin= ’ ,F11 . 6 , ’ ,  sf=  ’,fll.6,’  f rip= ’ ,f 11 . 6) 

write  (3,*)  ’next  is  total  pts.  and  no.  upstream  of  throat,’ 
write  (3,*)  ’then  x,y,  pe/pO,  dy/dx,  d2y/dx2  are:  ’ 
write  (3,*) 

>  ’waltan=dy/dx  and  sd=d2y/dx2  only  written  downstream  of  throat’ 

*  note  that  these  two  are  given  in  PERFC,  format  statement  89,  9-94  sps 

write  (3,*)  king+numpts,  numpts 

write  (3,103)  (xp(k) ,yp(k) ,pep0(k) ,k=l , numpts) 


write  (3,105)  (s(k) ,fs(k) ,spr(k) ,waltan(k) ,sd(k) ,k=l,king) 
103  f ormat (3(el4 . 7, lx) ) 

105  f ormat (5(el4 . 7, lx) ) 

else 

write  (3,*)  ’next  is  streamline  for  icall=  ’, icall 
write  (3,*)  ’totpts  and  no.  pts  upstream  of  throat,’ 
write  (3,*)  ’  then  x,y,  pe/pO:  ’ 

write  (3,*)  ’downstream  of  throat,  also  dy/dx  and  d2y/dx2’ 
write  (3,*)  king+numpts,  numpts 
write  (3,103)  (xp(k) ,yp(k) ,pep0(k) ,k=l , numpts) 
write  (3,105)  (s(k) ,fs(k) ,spr(k) ,waltan(k) ,sd(k) ,k=l,king) 
end  if 

* 

icall  =  icall  +  1 
return 

* 

*  following  all  format  statements  saved  for  comparison 

* 


*  write  (3,84)  RC , ETAD , AMACH , BMACH , CMACH , EMACH , MC , AH  476 

*  write  (3,90)  (K,S(K) ,FS(K) , WALT AN (K) ,SD(K) ,WMN(K) ,DMDX(K) ,SPR(K) ,D  478 

*  1PX(K) ,K=1,KING)  479 

C  493 

*  save  original  format  statements  from  perfc  for  reference  ******************** 

*84  FORMAT  (1H  ,4H  RC=,F11 . 6 , 3X , 5HETAD=F8 . 4 , 4H  DEG , 3X , 6HAMACH=F10 . 7 , 3X  494 

*  1 , 6HBMACH=F1 0 . 7 , 3X , 6HCMACH=F 1 0 . 7 , 3X , 6HEMACH=F 1 0 . 7 , 3X , A4 , 2HH=F 11.7/)  495 

*89  FORMAT  (1H  , 9X , 5HP0INT , 7X , 5HX (IN) , 9X , 5HY (IN) , 9X , 5HDY/DX , 8X , 7HD2Y/D  501 

*  1X2 , 7X , 8HMACH  NO . , 7X , 5HDM/DX , 9X , 5HPE/P0 , 1 IX , 6HDPR/DX/ )  502 

*90  FORMAT  (10 (10X , 13 , 2X , 0P6F14 . 7 , 1P2E16 . 5/) )  503 

*92  FORMAT  (1H  ,’  RC= ’ ,F11 . 7 , ’ ,  STREAMLINE  RATI0= ’, FI 1 . 8, ’ ,  TEST  505 

*  1  CONE  BEGINS  AT ’ , F12 . 7 , ’  IN . ’  /  )  506 

*95  FORMAT  (1H1,3A4,45H  INVISCID  NOZZLE  CONTOUR,  RADIAL  FLOW  ENDS  ATF1  509 

*  11.6, 25H  IN.,  TEST  CONE  BEGINS  ATF11.6.19H  IN.,  SCALE  FACT0R=F9 . 4/)  510 

END  527 


C  Sivells-to-Harris  Interface  Code 

*  PROGRAM  MAKEBLIN.F0R. 

*  Steven  P.  Schneider  Purdue  University  317-494-3343 

*  this  is  a  program  to  read  in  output  from  the  sivells  code, 

*  add  specifics  for  Re,  and  write  in  a  form  readable 

*  by  the  Harris  code  for  bl . 

*  specific  for  the  nozzle  block  problem  sps  6-90 

*  add  some  code  for  the  contraction  computation  12-5-90  sps 

*  add  code  for  output  of  file  for  arbitrary  shape  using  modified  Newtionian  thy 

*  sps  3-6-91 
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*  allow  for  arbitrary  power-law  distribution  of  points  3-6-91  sps, 

*  and  for  easy  change  of  j solve 

*  add  code  for  reading  in  pressure  data  from  euler  code,  and  shock  location 

*  data,  and  writing  an  output  file  for  the  harris  code  7-13-91  sps 

*  modified  to  fix  bug  with  error  code  on  different  case  7-28-91  sps 

*  modify  a  bit  to  use  with  whole  nozzle,  not  just  sivells  part  2-4-92  sps 

*  modify  to  make  print  stations  go  to  te  of  nozzle,  not  just  near.  4-8-93  sps 

*  modify  8-94  to  make  thermal  BC  explicit,  add  bullet  solution 

*  modify  9-8-94  to  skip  4  lines  in  sivells  for  new  sivells  radcur  format,  sps 

* 

parameter (maxpts=1000 , j  solve=l) 

*  j solve  is  the  number  of  solution  stations  per  station  of  input  data 

*  can  increase  to  make  finer  resolution  solution  possible 

character*80  text 

character *20  root , infile , outf ile ,parmf 1 , tawf ile 

*  most  of  these  are  the  number  of  stations  along  the  moc  soln 

*  ss  is  the  number  of  stations  in  bl  soln  to  get 

dimension  x(maxpts) ,y(maxpts) ,pepo (maxpts) ,s(maxpts) , theta (maxpts) 
dimension  proval (maxpts) ,prntval (maxpts) , ss ( j solve*maxpts) 
dimension  tw(maxpts) 

dimension  xsh (maxpts) ,ysh (maxpts)  ! shock  location 

common  /param/pstar  !for  passing  to  subroutines  for  arbitrary  shapes 
data  pi/3. 1415926535/, ksprnt/l/,ksprof/2/ 

*  * (print  info  every  ksprnt’th  soln  station;  print  profile  info 

*  *  (every  ksprof’th  soln  station) 

*  need  dense  printing  of  soln  for  gortler  work  to  get  good  values 

*  for  derivatives  of  wall  height  to  get  streamwise  curvature! 

*  need  lots  of  profiles  for  roughness  work  also 

* 

write  (*,*)  ’this  is  the  PC  version  of  makeblin’ 
write  (*,*)  ’unix  version  differs  in  namelist  format’ 
write  (*,*)  ’enter  root  filename  for  transfer:  ’ 
read  (*,5)  root 
5  format (a20) 

ileng  =  index (root,’  ’)  -  1 

write  (*,*)  ’enter  0  if  this  is  a  Sivells  (or  nozzle)  test,  ’ 

write  (*,*)  ’enter  1  if  this  is  a  flat  plate  test,  ’ 

write  (*,*)  ’enter  2  if  this  is  a  Lees  modified  newtonian  test:  ’ 

write  (*,*)  ’enter  3  if  Euler  output  for  body  is  to  be  read:  ’ 

write  (*,*)  ’enter  4  is  this  is  a  round  cone  at  zero  AQA:  ’ 

read  (*,*)  imodel 

if  (imodel  .eq.  0)  then 

write  (*,*)  ’Sivells  test  or  other  nozzle  test’ 
inf ile (1 : ileng)  =  root (1 : ileng) 
inf ile (ileng+1 : ileng+3)  =  ’.bl’ 

write  (*,*)  ’reading  input  data  from  file  ’, infile 
open(unit=l , f ile=inf ile , status=’ old’ ) 
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read  (1,10)  text 
10  format (a80) 

read  (1,10)  text 
read  (1,10)  text 
read  (1,10)  text 
read  (1,*)  numpts 

if  (numpts  .gt.  maxpts)  stop  ’too  many  points’ 
do  50  i=l, numpts 

read  (1,*)  x(i) ,y(i) ,pepo(i) 

50  continue 

close  (unit=l) 
write  (*,*) 

>  ’  wall.  temp.  BC?  Enter  1  if  isothermal,  2  adiabatic,  ’ 

write  (*,*)  ’  3  if  Tw/Taw  =  const:  ’ 
read  (*,*)  itflag 

write  (*,*)  ’read  itflag  as:  ’, itflag 

*  (done  reading  in  info  from  sivells  output  file) 

else  if  (imodel  .eq.  1)  then 

write  (*,*)  ’enter  numpts,  mach,  gam,  plend:  ’ 
read  (*,*)  numpts , amach, gam, plend 
gmexp  =  gam/ (gam-1) 
gmfact  =  (gam-l)/2.0 

denom  =  (1.0  +  gmf act*amach**2) **gmexp 

pepO  =  1.0/denom 

write  (*,*)  ’  gives  pep0=  ’ ,pep0 

do  60  i  =  1, numpts 

x(i)  =  plend*float(i-l)/float(numpts-l) 
y(i)  =1.0  !not  0.0,  messes  up  computations 
pepo(i)  =  pepO 
60  continue 

else  if  (imodel  .eq.  2)  then 

write  (*,*)  ’this  is  a  modified  newtonian  test;’ 
write  (*,*)  ’you  must  enter  the  shape  in  source  code’ 
write  (*,*)  ’enter  numpts,  mach:  ’  !,  pstar,  plend:  ’ 
read  (*,*)  numpts, amach  !, pstar, plend 
gam  =  1.4  lair 

if  (numpts  .gt.  maxpts)  stop  ’too  many  points’ 

*  now  compute  pressure  ahead  of  shock,  ratio  to  total  pressure 

*  in  stilling  chamber 

gmexp  =  gam/ (gam-1) 
gmfact  =  (gam-l)/2.0 

denom  =  (1.0  +  gmf act*amach**2) **gmexp 
pinfpO  =  1.0/denom 

write  (*,*)  ’  gives  pinf inity/p0=  ’, pinfpO 

*  now  compute  stagnation  or  total  pressure  behind  normal  shock, 

*  ratio  to  p_infty  ahead  of  shock(see  Anderson  p.  54,  3.17) 
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denoml  =  4*gam*amach**2  -  2* (gam-1) 

piecel  =  (  (gam+l)**2  *  amach**2/denoml  )** (gmexp) 

piece2  =  (l-gam+2*gam*amach**2)/(gam+l) 

ptpinf  =  piecel*piece2 

ptpO  =  ptpinf *pinfpO 

write  (*,*)  ’pstag  on  body/p  infty  is  1 , ptpinf 
write  (*,*)  ’gives  pstag  on  body/pstag  in  stilling= ’ ,ptpO 
cpmax  =  (2 . 0/ (gam*amach**2) )  *  (ptpinf  -  1)  ! Anderson  3.19 
write  (*,*)  ’cpmax  computed  as  ’, cpmax 

*  now  compute  body  shape 

*  for  non-sphere  shape,  do  this  in  subroutine 

c  write  (*,*)  ’passing  to  blunts;  pstar ,plend=’ ,pstar ,plend 

*  blunts  commented  out  12-93,  find  subroutine  to  bring  back 

c  call  blunts(plend,numpts,x,y, theta) 

*  following  does  a  sphere 

do  70  i  =  l.numpts 

x(i)  =  -1.0  +  float (i-1) /float (numpts-1) 
y (i)  =  sqrt (1.0  -  (x(i))**2) 
if  (y(i)  .gt.  0.0)  then 

dydx  =  -l*x(i)/y(i)  (needed  for  newtonian  thy 
theta(i)  =  atan(dydx) 
else 

theta(i)  =  pi/2.0 
end  if 
70  continue 

write  (*,*)  ’last  x,y  are  ’ ,x(numpts) , ’ , ’ ,y(numpts) 
do  80  i  =  l.numpts 

*  now  compute  ratio  of  pe  to  ptotal  ahead  of  shock,  using  Lees 

*  modified  newtonian  thy  -  formula  derived  using  p.  54 

pepo(i)  =  (ptpO  -  pinfpO) * (sin(theta(i) ) ) **2  +  pinfpO 
if  (pepo(i)  .ge.  1.0)  pepo(i)  =  0.99999  Iso  not  singular 
80  continue 

else  if  (imodel  .eq.  3)  then 

write  (*,*)  ’working  for  euler  data  for  blunt  body’ 
write  (*,*)  ’enter  mach,  gamma:  ’ 
read  (*,*)  amach,gam 

*  now  compute  pressure  ahead  of  shock,  ratio  to  total  pressure 

*  in  stilling  chamber 

gmexp  =  gam/ (gam-1) 
gmfact  =  (gam-l)/2.0 

denom  =  (1.0  +  gmf act*amach**2) **gmexp 
pinfpO  =  1.0/denom 

write  (*,*)  ’  gives  pinf inity/p0=  ’, pinfpO 

*  now  compute  stagnation  or  total  pressure  behind  normal  shock, 

*  ratio  to  p_infty  ahead  of  shock(see  Anderson  p.  54,  3.17) 

denoml  =  4*gam*amach**2  -  2* (gam-1) 

piecel  =  (  (gam+l)**2  *  amach** 2 /denoml  )** (gmexp) 
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piece2  =  (l-gam+2*gam*amach**2)/(gam+l) 
ptpinf  =  piecel*piece2 
ptpO  =  ptpinf *pinfpO 

write  (*,*)  ’pstag  on  body/p  infty  is  1 , ptpinf 

write  (*,*)  ’gives  pstag  on  body/pstag  in  stilling= ’ ,ptpO 

write  (*,*)  ’reading  euler  data  for  blunt  body:  ’ 

inf ile (1 : ileng)  =  root (1 : ileng) 

inf ile (ileng+1 : ileng+4)  =  ’.eul’ 

write  (*,*)  ’reading  input  data  from  file  ’, infile 
open(unit=l , f ile=inf ile , status=’ old’ ) 
read  (1,90)  lineskip 

90  format (i6)  Inumber  of  lines  to  skip 

do  92  i  =  1, lineskip  (skip  ’lineskip’  lines  of  text 
read  (1,91)  text 

91  format (a80) 

92  continue 

read  (1,*)  numpts 

if  (numpts  .gt.  maxpts)  stop  ’too  many  points’ 
do  95  i=l, numpts 

read  (1,*)  x(i) ,y (i) ,pepinf  Iread  euler  data  for  pressure 
pepo(i)  =  pepinf *pinfp0 

95  continue 

*  Now  read  shock  location  data: 
read  (1,90)  lineskip 

do  96  i  =  1, lineskip  ! skip  ’lineskip’  lines  of  text 
read  (1,91)  text 

96  continue 

read  (1,*)  numshpts 

if  (numshpts  .gt.  maxpts)  stop  ’too  many  points’ 
do  97  i=l, numshpts 

read  (1,*)  xsh(i) ,ysh(i)  Iread  euler  data  for  shock  location 

97  continue 

write  (*,*)  ’done  reading  from  file’ 

else  if  (imodel  .eq.  4)  then 

write  (*,*)  ’  cone:  enter  numpts,  gam,  axial  length:  ’ 
read  (*,*)  numpts ,gam, axleng 

write  (*,*)  ’  enter  half -angle  (deg.),  shock  angle,  pepO:  ’ 
read  (*,*)  anghalf ,  wave,  pepO 
anghalf  =  anghalf *pi/180 . 0 
do  100  i  =  1, numpts 

x(i)  =  axleng*float (i-1) /float (numpts-1) 

y(i)  =  tan (anghalf )*x(i)  Inot  0.0,  messes  up  computations 

pepo(i)  =  pepO 
100  continue 
else 

stop  ’imodel  must  be  0,1, 2, 3,  or  4’ 
end  if 
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*  now  check  computed  values 

do  150  i  =  l,numpts 

if  (imodel  .eq.  2  .or.  imodel  .eq.  3)  then 

if  (pepo (i) -ptpO  .gt.  -0.01*ptp0)  then  ! arbitrary  nearness 
write  (*,*)  ’pepo( 1 , i , 1 )=  ’,pepo(i) 
write  (*,*)  ’ptp0=  ’ ,ptp0 

write  (*,*)  ’if  first  gt  second  when  computed  by  VGBLP,’ 
write  (*,*)  ’this  will  give  surface  static  pressure  larger 
write  (*,*)  ’  than  the  total  pressure  on  the  surface’ 
write  (*,*)  ’this  would  be  a  fatal  error  in  VGBLP’ 
if  (pepo(i)  .gt.  ptpO)  then 

write  (*,*)  ’reducing  pepo(’,i,’)  to  0 . 99999*ptp0 ’ 
write  (*,*)  ’  to  forstall  error  in  VGBLP’ 
pepo(i)  =  0.99999*ptp0 
end  if 
end  if 
end  if 
150  continue 

* 

*  now  read  in  parametric  info  from  parameter  file 

* 

parmf 1 (1 : ileng)  =  root (1 : ileng) 
parmfl(ileng+l : ileng+3)  =  ’.re’ 

write  (*,*)  ’reading  reynolds  number  scaling  info  from  ’, parmf 1 

open(unit=l , f ile=parmf 1 , status= ’ old’ ) 

read  (1,5)  outfile  !read  filename  to  write  to 

read  (1,*)  throat  ! throat  radius  in  feet 

*  (assumes  input  scaled  so  throat  radius  is  one  unit)  * 

*  (for  more  general  shapes,  just  treats  ’throat’  as  scaling  parameter 

*  for  lengths) 

read  (1,*)  prandtl 

*  ptotal,ttotal,xmach  are  conditions  at  infinity  -  see  p.  44 

*  for  nozzle  are  stagnation  chamber  conditions 

*  for  non-nozzle,  xmach  seems  to  affect  mostly  the  computations 

*  involving  the  flow  behind  the  shock.  Should  be  freestream  values! 

read  (1,*)  ptotal 
read  (1,*)  rgas 

read  (1,*)  ttotal  ! ahead  of  le  shock 
read  (1,*)  xmachi 
if  (imodel  .ne.  0)  then 
if  (xmachi  .It.  1.0)  then 

write  (*,*)  ’xmach  given  as  ’, xmachi 

write  (*,*)  ’should  be  freestream  value  ahead  of  shock,’ 
write  (*,*)  ’  not  the  value  at  stagnation!!’ 
end  if 
end  if 

close  (unit=l) 
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* 

*  now  set  defaults  for  input  to  harris  program  (besides  harris’s) 

* 

if  (imodel  .eq.  2  .or.  imodel  .eq.  3)  then 
ibody  =  1  ! stag  pt  at  nose 

j  =  1  ! axisymmetric 

else  if  (imodel  .eq.  1)  then 

ibody  =2  !no  stagnation  point  at  nose 
j  =  0  !  2D 

write  (*,*)  ’assuming  2D  geometry1 
else  if  (imodel  .eq.  0)  then 

write  (*,*)  ’enter  0  if  2D  nozzle,  1  if  axisymmetric:  ’ 
read  (*,*)  j 

if  (j  .ne.  0  .and.  j  .ne.  1)  stop  ’j  must  be  0  or  1’ 
else  if  (imodel  .eq.  4)  then 

ibody  =  2  (sharp  cone,  no  stag,  pt . 

j  =  1 

end  if 

ie  =  51  (from  test  case  number  4 

if  (imodel  .eq.  3)  then 

ientro  =  2  (variable  entropy  calculation 
else 

ientro  =  1 
end  if 

igeom  =  1  (create  coords  using  geometric  series;  need  xk,ie,xend(! 
kodunit=0  (US  units 

if  (itflag  .eq.  1  .or.  imodel  .eq.  4)  then 

kodwal=l  (specify  wall  temperature  distribution  (no  time  to  heat) 

write  (*,*)  ’  setting  kodwal=l,  isothermal  wall!!’ 

write  (*,*)  ’setting  wall  temp  equal  to  total  temp,  beware!’ 

*  model  temp  same  as  total  temp,  ambient 

*  note  that  total  temperature  behind  shock  is  the  same  as  ahead  of  shock 

*  *(not  quite  stagnation  point  temp,  but  close)* 
twall  =  ttotal 

else  if  (imodel  .eq.  2)  then 
kodwal=l 

write  (*,*)  ’enter  isothermal  wall  temperature,  rankine :  ’ 
read  (*,*)  twall 

write  (*,*)  ’using  wall  temp=  ’, twall 
else  if  (itflag  .eq.  2  .or.  imodel  .eq.  3)  then 
kodwal=2  (specify  adiabatic  wall 

write  (*,*)  ’  setting  kodwal=2,  adiabatic  wall!!’ 
else  if  (itflag  .eq.  3)  then 
kodwal  =  1 

write  (*,*)  ’setting  kodwal=l,  isothermal  wall’ 

write  (*,*)  ’enter  file  to  read,  .prt  file  with  Taw/TTl  data:  ’ 

read  (*,5)  tawfile 

write  (*,*)  ’reading  Taw/TTl  from  ’, tawfile 
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write  (*,*)  ’enter  const,  where  Tw/Taw  =  const:  ’ 
read  (*,*)  tratio 

write  (*,*)  ’using  tratio=  ’, tratio 
open  (unit=ll , f ile=tawf ile , status=’ old’ ) 
read  (11,10)  text 
read  (11,10)  text 

read  (11,*)  numptsl .prandtll ,j 1 , omegal ,rrefl ,urefl 
if  (numpts  .ne.  numptsl)  then 

numpts  =  numptsl  Ireset  and  solve  at  prev.  solution  stations! 
iendl  =  numpts 

throat  =  1.0  Idon’t  rescale  dimensions!! 
else 

iendl  =  j solve*numpts 
end  if 

do  170  i  =  1, numpts 

read  (11,*)  zl,rmil,sl,yel,dltastl,thetal,resl,pel,tel,uel, 

>  twott 1 , amache 1 , amue 1 , xi 1 , qsdl , hdl 

tw(i)  =  ttotal*twottl*tratio 
x(i)  =  zl 
y(i)  =  rmil 
pepo(i)  =  pel/ptotal 
170  continue 

close  (unit=ll) 
else 

stop  ’logic  error  in  setting  wall  thermal  BC’ 
end  if 

if  (itflag  .ne.  3)  then  ! otherwise  set  above  when  read  file 
iendl  =  jsolve*numpts  ! number  of  soln  stations 
end  if 

proinc  =  10.0  !  hopefully,  none 
prntinc  =  10.0 

sst  =  le+20  !no  transition  on  body  until  then  (laminar  flow) 
if  (imodel  .gt.  1  .and.  imodel  .It.  4)  then 

wave  =  90.0  ! shock  wave  angle  at  s=0,  needed  for  test  case  4  type  flows 

else  if  (imodel  .eq.  4)  then 

write  (*,*)  ’read  shock  wave  angle  at  origin  as:  ’ ,wave 
else 

wave  =  0.0  ! needed  for  shockless  type  flows 

end  if 

c  xend  =  120  !from  blasius  test  case 
xend  =10  ! as  in  test  case  4 

c  xk  =  1.275  lvalue  used  in  test  cases  in  book,  sets  grid 

c  xk  =  1 . 1  ! because  value  used  in  text  gives  hyper-dense 

*  grid  near  wall,  which  makes  for  difficulties, 

c  xk  =  1.0  Hike  test  case  number  4 

c  xk  =  1.05  ! because  1.0  gives  little  near  wall  for  sphere 

xk  =  1.1  ! because  1.05  gives  not  great  resolution  for  stability 

* 
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*  compute  the  arc  length  along  the  wall  (see  (66)  of  harris  paper) 

*  (approximate  with  straight  line  segments  between  stations) 

s ( 1 )  =  0.0 

if  (itflag  .le.  2)  then 
tw(l)  =  twall 
end  if 

do  180  i  =  2,numpts 

s (i)  =  s(i-l)  +  sqrt (  (x(i)-x(i-l))**2  +  (y (i) -y (i-1) ) **2  ) 
if  (itflag  .le.  2)  tw(i)  =  twall  ! set  above,  assumes  same  along  wall 
180  continue 

*  compute  the  s  stations  to  get  soln  at,  and  to  write  at,  even  spacing 

*  in  sqrt(s)  normally,  other  power  for  other  cases 

if  (imodel  .It.  2)  then 
power  =  2.0 

write  (*,*)  3 using  square  root  distribution  of  pts3 
else 

*  with  sphere,  have  problems  with  stepsize  increasing  too  rapidly  near  le 

power  =  1.0 

write  (*,*)  3  using  linear  distribution  of  pts,  good  for  blunt 3 

*  *(try  this,  since  problems  with  T<0  at  le)*like  test  case  4 
end  if 

rootl  =  1.0/power 

rsinc  =  (s (numpts) -s (1) ) **rootl/f loat (iendl)  !try- 

ss(l)  =  rsinc**power 

ss(2)  =  ss(l)  (required 

ss(3)  =  ss(l)  !also  required,  actually 

srun  =  3*ss(l)  (running  value  of  s 

iprnt  =  0 

ipro  =  0 

do  200  i  =  4, iendl 
srunold  =  srun 

srun  =  (  (srunold) **rootl  +  rsinc  )**power 
ss(i)  =  srun-srunold 

*  add  last  part  to  following  to  have  a  print  station  at  the  end: 

if  (mod(i ,ksprnt)  .eq.  0  .or.  i  .eq.  1  .or.  i  .eq.  iendl)  then 
iprnt  =  iprnt  +  1 
if  (iprnt  .gt.  maxpts)  then 

write  (*,*)  3 iprnt  =  iprnt, 3  exceeded  maxpts3 
stop  ’fatal  error3 
end  if 

prntval (iprnt)  =  srun 
end  if 

if  (mod(i ,ksprof )  .eq.  0)  then 
ipro  =  ipro  +  1 
proval(ipro)  =  srun 
end  if 
200  continue 

* 
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*  Now  write  out  NAM1  namelist  into  file 

* 

write  (*,*)  ’  writing  bl  input  to  file  ’.outfile 

open(unit=2 , f ile=outf ile) 

write  (2,*)  ’&NAM1 1 

write  (2,*)  ’ IBQDY= ’ , ibody 

write  (2,*)  ’IE=’,ie 

write  (2,*)  ’ IEND1= ’ , iendl 

write  (*,*)  ’may  need  to  reset  param.inc,  note  that  ’ 

write  (*,*)  ’ie  is  now  ’,ie,’  and  iendl=  iendl 

write  (2,*)  ’ IENTRD=’ , ientro 

write  (2,*)  ’ IGE0M= ’ , igeom 

write  (2,*)  ’IPR0=’,ipro 

write  (2,*)  ’ IPRNT= ’ , iprnt 

write  (*,*)  ’ ipro  and  iprnt  are  ’,ipro, iprnt 

write  (2,*)  ’J=’,j 

write  (2,*)  ’KDDUNIT=’ ,kodunit 

write  (2,*)  ’KDDWAL=’ ,kodwal 

if  (imodel  .ne.  2)  then 

phii  =  (180 . 0/pi) *atan(  (y(2)-y(l))/(x(2)-x(l))  ) 
else 

write  (*,*)  ’setting  leading  edge  angle  =  90  degrees’ 
phii  =90.0 


end  if 

write 

(2,*) 

’ PHII= ’ ,phii 

write 

(2,*) 

’PR=’ .prandtl 

write 

(2,*) 

’ PRNTINC= ’ .prntinc 

write 

(2,*) 

’PRNTVAL=’ , (prntval(i)*throat,i=l, iprnt) 

write 

(2,*) 

’PR0INC=’ .proinc 

write 

(2,*) 

’PR0VAL=’ , (proval (i) *throat , i=l , ipro) 

write 

(2,*) 

’PT1=’ ,ptotal 

write 

(2,*) 

’ R= ’ , rgas 

write 

(2,*) 

’SST=’ ,sst 

write 

(2,*) 

’TT1=’ ,ttotal 

write 

(2,*) 

’ WA¥E= ’ .wave 

write 

(2,*) 

’ XEND= ’ , xend 

write 

(2,*) 

’ XK= ’ ,xk 

write 

(2,*) 

’ XMA= ’ , xmachi 

write 

(2,*) 

’/’  lend  of  namelist  input 

*  end  of  writing  : 

naml .  Now  compute  and  write  nam2: 

write 

(2,*) 

’  &NAM2  ’ 

write 

(2,*) 

’NUMBERS .numpts 

write 

(*,*) 

’and  NUMBER  is  ’.numpts 

write 

(2,*) 

’PE=’ , (pepo(i)*ptotal,i=l, numpts) 

write 

(2,*) 

’RMI=’ , (y (i) *throat , i=l .numpts) 

write 

(2,*) 

’ S= ’ , (s (i) *throat , i=l .numpts) 

write 

(2,*) 

’ SS=’ , (ss (i) *throat , i=l , iendl) 

if  (kodwal 

.eq.  1)  then  ! specified  wall  temp 
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write  (2,*)  ’TW=’ , (tw(i) , i=l ,numpts) 
else  ! specified  heat  transfer 

write  (2,500)  numpts-1  Imake  adiabatic  wall 

500  format (’QW=’ ,i4, ’*0.0,0. O’)  ! , 0 . 0  kluge  for  compiler  bug 
end  if 

write  (2,*)  ’ Z= ’ , (x(i) *throat , i=l ,numpts) 
write  (2,*)  ’/’ 

*  end  of  writing  nam2.  Now  write  nam3,  if  required 

if  (imodel  .eq.  3)  then  ! doing  entropy  computations  through  shock 
write  (2 , *)  ’&NAM3’ 
write  (2,*)  ’NUMBERS,  numshpts 
write  (*,*)  ’now  numshpts  (JL)  is  1 , numshpts 
write  (2,*)  ’RRS=  ’ , (xsh(i) *throat , i=l , numshpts) 
write  (2,*)  ’ZZS=  ’ , (ysh(i) *throat , i=l , numshpts) 
end  if 

write  (2,*)  ’/’ 
close  (unit=2) 

* 

stop  ’end  of  makeblin’ 
end 

D  Harris  to  e**MALIK  Translation  Code 

*  BLT0STAB.F0R 

*  this  is  a  program  to  take  vgblp  data  and  put  in  E**MALIK  form 

*  sps  7-2-90 

*  revised  7-9-90  to  fit  with  input  form  provided  by  Maliks  vgblp  program, 

*  different  from  that  implied  in  preliminary  paper 

*  revised  7-17-90  to  fix  problem  with  getting  correct  matches 

*  also  revised  to  give  correct  scaling  of  profiles  for  emalik  code 

*  revised  7-18-90  to  scale  before  passing  to  utder,  makes  cutoffs  clearer 

*  revised  7-16-91  to  write  more  info  to  .cur  file  sps 

*  revised  7-16  and  7-17-91  to  get  better  derivatives  of  surface  shape  sps 

*  and  also  to  skip  past  unused  iterations  printed  using  variable  entropy 

*  computations 

*  revised  8-27-91  to  use  derivative  data  now  output  by  Harris  code 

*  harris  code  outputs  FZ  and  TZ,  which  is  actually  the  derivs  wrt  the 

*  wall-normal  coordinate  eta  at  the  next  streamwise  solution  station. 

*  However,  we  do  not  normally  output  every  single  solution  station,  so 

*  cannot  just  use  this  at  the  next  solution  station  output.  So  accept 

*  the  error  involved  in  using  the  derivs  at  next  solution  station  in 

*  place  of  derivs  at  current,  for  now. 

*  mod.  5-27-94  to  read  and  compute  derivs  in  double,  needed  for  t  profiles 

*  mod.  9-2-94  sps  to  read  derivs  of  nozzle  contour  from  .bl  sivells  file 

*  (couldnt  get  good  derivatives  of  nozzle  contour  from  printed  data) 

*  mod  9-9-94  to  interpolate  sivells  data,  stations  don’t  match,  sps 
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*  mod.  10-14-94  to  read  curvature  data  from  file  when  not  using  sivells,  sps 

*  mod  10-24-94  to  use  Frank  Chen’s  method  of  getting  curvature  from  nozzle 

*  contour  data,  using  old  NOS  routine,  sps. 

* 

c  implicit  double  precision  (a-h,o-z)  Ineeded  for  derivs 
parameter  (mpt=1000 ,npt=1000) 

character*81  text  !try  changing  to  81  from  80  per  7-91  version 
character *20  blf ile , outf il ,root ,prof il , curf il , sivf il ,ref il ,gortf 1 
character*20  sivck,unsivck 

dimension  z(mpt) ,rmi(mpt) ,s(mpt) ,ye(mpt) ,dltast(mpt) ,theta(mpt) , 

>  res(mpt) ,pe(mpt) ,te(mpt) ,ue(mpt) ,twttl(mpt) ,ame(mpt) , 

>  amue (mpt ) , xi (mpt ) 

dimension  xsiv(mpt) ,ysiv(mpt) ,ssiv(mpt) ,dydx(mpt) , 

>  d2ydx2(mpt) ,radcur(mpt) 

dimension  df (mpt) , coef (mpt ,4) ,wk(7*mpt+9) ,rrmi (mpt)  !for  CSDS  subroutine 
dimension  x(npt) ,u(npt) ,ul(npt) ,u2(npt) ,t(npt) ,tl(npt) ,t2(npt) 
dimension  uldeb(npt) , etascal (npt) ,t23b(npt) 
data  small/1 . 0e-4/ 

data  norder/3/  ! order  of  polynomial  interpolation 

data  maxiter/5/  (maximum  number  of  variable  entropy  iterations  in  file 
data  numiter/0/  Inumber  of  variable  entropy  iter  written  to  file 
data  nhophill/0/  lused  as  an  offset 

data  eps/0.1/  ! fractional  error  acceptable  in  derivs,  this  is 

*  really  a  check  on  the  unit  conversions 
data  huge/1. OelO/  ! if  radcur  infinite 

data  df /mpt*l . e-3/  (estimate  of  standard  dev.  for  CSDS 

* 

write  (*,*)  ’enter  root  filename  for  writing  and  reading:  ’ 
read  (*,5)  root 
5  format (a20) 

il  =  index (root,’  ’)  -1 
curfil(l:il)  =  root(l:il) 
curf il(il+l : il+4)  =  ’.cur’ 

write  (*,*)  ’opening  ’.curfil,’  for  curvature  output’ 
open(unit=3 , f ile=curf il) 

write  (3,*)  ’curvature  data  for  this  case’ 
write  (3,*)  ’first  set:  nz,z,s,rmi, radcur, epsxr’ 
gortfl(l:il)  =  root(l:il) 
gortf 1 (il+1 : il+4)  =  ’.gor’ 

write  (*,*)  ’opening  ’.gortfl,’  for  gortler  number  output’ 
open(unit=7 , f ile=gortf 1) 

write  (7,*)  ’z,s, theta, res, radcur, gortno  for  file:  ’ ,root 

sivck(l:il)  =  root(l:il) 

sivck(il+l : il+4)  =  ’.sck’ 

unsivck(l : il)  =  root(l:il) 

unsivck(il+l : il+4)  =  ’.usk’ 

blfile(l:il)  =  root(l:il) 

blf ile(il+l : il+4)  =  ’.prt’ 
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write  (*,*)  ’  opening  ’, bifile,’  for  read  printed  info’ 
open(unit=l , f ile=blf ile , status= ’ old1 ) 
read  (1,10)  text 
10  format (a80) 

read  (1,10)  text  ! second  line  of  text 

*  for  files  with  variable  entropy  computations,  several  iterations 

*  may  exist  in  the  file,  so  the  print  data  is  redone  for  several 

*  iterations,  following  sequential  in  file.  Skip  to  last 

*  read  iprnt ,  prandtl,  j  which  stays  the  same 

read  (1,*)  iprnt, prandtl, jgeom, omega, rref.uref 
if  (iprnt  .gt.  mpt)  stop  ’too  many  print  stations’ 

*  now  read  the  first  set  of  iprnt  values: 

do  110  iiter  =  l,maxiter 
do  100  i  =  1 , iprnt 

read  (1,*)  z(i) ,rmi(i) ,s(i) ,ye(i) ,dltast(i) ,theta(i) ,res(i) , 

>  pe(i) ,te(i) ,ue(i) ,twttl(i) ,ame(i) ,amue(i) ,xi(i) 

if  (z(i)  .It.  0.0)  nupstrm  =  i  Inumber  of  pts  upstream  of  throat 
100  continue 

read  (1 , * , end=120)  iprnt2 ,prandtl2 , jgeom2 , omega2  Inow  repeated 
read  (1 , * , end=120)  z(l) ,rmi(l) ,s(l) , 

>  ye(l) ,dltast(l) ,theta(l) ,res(l) , 

>  pe(l),te(l),ue(l) ,twttl(l) ,ame(l) ,amue(l) ,xi(l) 

if  (s(l)  .It.  s (iprnt))  then  ! there  IS  another  iteration 
numiter  =  numiter  +  1 

write  (*,*)  ’read  iteration  ’, numiter,’  read  next’ 
backspaced)  Iback  up  to  before  the  first  line  in  this  iter, 
else 

write  (*,*)  ’iprnt  =  ’, iprnt 

write  (*,*)  ’s(l)=  ’,s(l),’  s(iprnt)=  ’,s(iprnt) 
stop  ’fatal  error,  funny  business  in  reading  print  file’ 
end  if 
110  continue 

stop  ’reached  maxiter  reading  variable  entropy  data’ 

120  continue  ! reached  eof  looking  for  next  iteration,  done- 
close  (unit=l) 

write  (*,*)  ’there  are  ’ ,numiter+l ,’ iteration  sets  in  file’ 

*  now  open  profile  info  file: 

profil(l:il)  =  root(l:il) 
prof il (il+1 : il+4)  =  ’.pro’ 

write  (*,*)  ’opening  ’,profil,’  for  reading  from  vgblp’ 
open(unit=l , f ile=prof il , status= ’ old’ ) 
read  (1,150)  text 
150  format (a80) 

read  (1,*)  ipro 

* 

outfil(l:il)  =  root(l:il) 
outfil (il+1 : il+4)  =  ’.bfl’ 

write  (*,*)  ’opening  ’, outfil,’  for  writing  to  e**malik’ 
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open(unit=2 , f ile=outf il , form=’ unformatted’ ) 
write  (2)  text 
write  (2)  ipro 

* 

*  Decide  if  getting  curvature  data  from  sivells  or  elsewhere: 

write  (*,*)  ’where  to  get  curvature  data  for  gortler?  ’ 
write  (*,*)  ’enter  1  if  sivells  data,  2  if  by  spline  diff . :  ’ 
read  (*,*)  iwherec 

* 

*  beginning  of  block  where  use  1  of  2  ways  to  get  curvature  data 

* 

if  (iwherec  .eq.  1)  then 

*  now  open  the  .re  file  to  get  scaling  information  for  the  sivells  file 

refil(l:il)  =  root(l:il) 
ref il (il+1 : il+3)  =  ’.re’ 
write  (*,*)  ’opening  ’,refil, 

>  ’  to  get  reynolds  scaling  for  sivells’ 

open  (unit=4 , f ile=ref il , status= ’ old ’ ) 
read  (4,150)  text  ! skip  first  line 
read  (4,*)  throatrad  ! in  feet,  scales  sivells  data 
close  (unit=4) 

* 

*  now  open  the  sivells  output  file  directly,  the  XXXX.bl  file  used 

*  as  input  to  the  harris  code.  Pick  up  the  derivatives  of  the  nozzle 

*  contour  from  here.  9-94  sps 

* 

sivfil(l:il)  =  root(l:il) 
sivfil (il+1 : il+3)  =  ’.bl’ 

write  (*,*)  ’opening  ’, sivfil,’  to  read  sivells  contour  derivs’ 

open  (unit=4,f ile=sivf il , status=’ old’ ) 

read  (4,150)  text 

read  (4,150)  text 

read  (4,150)  text 

read  (4,150)  text  ! skip  header  lines 

read  (4,*)  ntotal,nhophill  ! total  num  pts,  no.  of  hopkins-hill 
do  160  i  =  l,nhophill 

read  (4,*)  xdum , ydum , pratdum  ! skip  past  these  points 
if  (i  .eq.  1)  then 
sarcl  =  0.0 
xold  =  xdum 
yold  =  ydum 
else 

sarcl  =  sarcl  +  sqrt ( (xdum-xold) **2  +  (ydum-yold) **2) 
xold  =  xdum 
yold  =  ydum 
end  if 
160  continue 

nsiv  =  ntotal-nhophill 
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do  200  i  =  l,nsiv 

read  (4,*)  xsiv(i) ,ysiv(i) ,pratdum,dydx(i) ,d2ydx2(i) 
xsiv(i)  =  xsiv(i)*throatrad 
ysiv(i)  =  ysiv(i)*throatrad 
if  (i  .eq.  1)  then 

ssiv(i)  =  throatrad*sarcl 
else 

ssiv(i)  =  ssiv(i-l)  + 

>  sqrt ( (xsiv(i)-xsiv(i-l) )**2  +  (ysiv(i)-ysiv(i-l))**2) 

end  if 

d2ydx2(i)  =  d2ydx2(i)/throatrad  (change  to  feet  from  throat  radii 
if  (d2ydx2(i)  .ne.  0.0)  then 

radcur(i)  =  (1 . 0+dydx(i) **2) **1 . 5/d2ydx2 (i) 
else 

radcur(i)  =  huge 
end  if 

*  (concave  is  minus  for  malik! 

200  continue 

close  (unit=4) 

*  (note  that  must  be  checked  that  these  are  at  same  stations) 

write  (*,*)  ’opening  ’,sivck,’  for  sivells  curvature  check’ 
open  (unit=4,f ile=sivck) 

write  (4,*)  ’sivells  curvature  check:  i,xsiv,ysiv,ssiv,radcur’ 
do  210  i  =  l,nsiv 

write  (4,209)  i,xsiv(i) ,ysiv(i) ,ssiv(i) ,radcur(i) 

209  f ormat (i4,3(lx,f 12 .5) , lx, lp , el2 . 5) 

210  continue 
close  (unit=4) 

ngor  =  iprnt  -  nupstrm  (number  of  pts  in  .prt  file  downstream  of  throat 

* 

else  if  (iwherec  .eq.  2)  then  (get  derivs  using  contour  directly 

* 

write  (*,*)  ’opening  ’ ,unsivck, ’  for  unsivells  curv.  check’ 
open  (unit=4,f ile=unsivck) 

write  (4,*)  ’non-sivells  curv.  ck:  z , rrmi , dydx , d2ydx2 , radcur ’ 
c  From  Frank_Chen.AERONAUTICS@qmgate.larc.nasa.gov  Mon  Oct  17  13:53  EST  1994 
c  RE>gortler  test  case.  The  fragment  of  the  code  I  used  is  very  simple. 

*  following  uses  NOS  routine  CSDS,  see  header  for  this  subroutine. 

IPT1=-1 

fnmin  =  iprnt  -  (2 . 0*iprnt) **0 . 5 

fnmax  =  iprnt  +  (2 . 0*iprnt) **0 . 5 

fn  =  (fnmin+fnmax) *0 . 5  (a  guess  for  what  to  use 

CALL  CSDS (mpt , iprnt , Z , RMI , DF , f n , IPT1 , C0EF , WK , IERR) 

if  (ierr  .ne.  0)  then 

write  (*,*)  ’error  return  from  CSDS,  ierr=  ’,ierr 
stop  ’halting’ 
end  if 

rrmi(l)  =  coef(i,l)  ! dh  =  0  for  these  three,  the  first  point 
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if  (rrmi(l)  .eq.  0.0)  then 

write  (*,*)  ’problem  with  csds  at  first  pt . ,  rrmi(l)=0’ 
write  (*,*)  ’set  to  rmi(l)’ 
rrmi ( 1 )  =  rmi ( 1 ) 
end  if 

dydx(l)  =  coef(i,2) 
d2ydx2(l)  =  2 . 0*coef (i , 3) 

DO  217  1=1 , iprnt-1 
DH=Z (1+1) -z (i) 

*  *  note  that  rrmi  is  the  interpolated  location  of  rmi  from  spline  fit* 
RRMI ( 1+1 ) = ( (COEF (1,4) *DH+C0EF (1,3)) *DH+C0EF (1,2)) *DH+C0EF (1,1) 
dydx(I+l)=(3 . 0*C0EF(I ,4) *DH+2 . 0*C0EF(I , 3) ) *DH+C0EF(I , 2) 

d2ydx2 (1+1) =6 . 0*C0EF (I ,4) *DH+2 . 0*C0EF (1,3) 

217  CONTINUE 

*  end  of  frank  chen  fragment  (which  has  been  adapted  here) 

do  220  i  =  1 , iprnt 

if  (d2ydx2(i)  .ne.  0.0)  then 

radcur(i)  =  (1 . 0+dydx(i) **2) **1 . 5/d2ydx2 (i) 
else 

radcur(i)  =  huge 
end  if 

write  (4,219)  z(i) ,rrmi(i) ,dydx(i) ,d2ydx2(i) ,radcur(i) 

219  format (5(lx, lp,el4 . 7) ) 

*  ! concave  is  minus  for  malik! 

220  continue 
close  (unit=4) 

ngor  =  iprnt  !for  gortler  printout 
else 

stop  ’invalid  iwherec’ 
end  if 

* 

*  end  of  block  where  get  radcur  in  one  of  two  ways 

* 

*  now  write  gortler  number  output  for  checking 
do  225  i  =  l,ngor 

if  (iwherec  .eq.  2)  then  !s  array  and  radcur  array  indexed  same 
retheta  =  res(i)*theta(i)/s(i) 
radcur 1  =  radcur (i) 
i2  =  i 

else  ! using  sivells  output 
i2  =  i  +  nupstrm 

si  =  s(i2)  !i  indexes  over  .prt  array,  NOT  sivells  array 

call  locate(ssiv,nsiv,sl , jsiv) 
if  (jsiv  .ge.  norder)  then 

call  polint (ssiv(jsiv-norder+l) , radcur (jsiv-norder+1) , 

>  norder, si, radcurl , errest) 

else 

call  polint (ssiv(l) , radcur (1) , norder , 
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>  sl,radcurl,errest) 
end  if 

if  (s(i2)  .eq.  0.0)  then 

write  (*,*)  ’i2,s=  ’ ,i2,s(i2),’  s(i2-l)=  ’,s(i2-l) 
stop  ’fatal’ 
end  if 

retheta  =  res (i2) *theta(i2) /s (i2) 
end  if 

if  (radcurl  .It.  0.0)  then 

gortno  =  retheta  *  sqrt(theta(i2)/abs(radcurl)) 
else  if  (radcurl  .gt.  0.0)  then 
gortno  =  0.0  ! convex 

else 

gortno  =  huge 
end  if 

write  (7,223)  z(i2) ,s(i2) ,theta(i2) ,res(i2) .radcurl, gortno 
223  format (6(lx, lp, el2 . 5) ) 

225  continue 
write  (7,*) 

>  ’now  z(in) ,s(ft) ,radcur,gortl, gortno  from  emalik  algorithm:  ’ 

* 

*  Now  have  everything  need  from  prntval  stations.  Start  reading 

*  data  from  proval  stations  and  writing  to  malik  program 

* 

*  Before  read  the  profiles,  skip  to  the  last  set  of  profiles 

*  (multiple  sets  if  doing  variable  entropy  computations) 

*  Know  how  many  iterations  in  file  from  prt  file,  use  this  info  here 
do  250  iskip  =  l,numiter  Inumiter  is  num  in  file  -1 

do  240  istation  =  1 , ipro  ! ipro  stations 
read  (1,*)  nnp,sl 
do  230  j  =  1 , nnp 

read  (1,*)  xdum,udum,tdum,ulndum,tlndum 
continue 
continue 
continue 

open(unit=4, f ile=’bltostab . deb’ ) 

do  1000  nz  =  1 , ipro  Hoop  over  stations 

first,  get  general  info  for  station  from  iprnt  file 

read  (1,*)  nnp, si  Inumber  of  points  in  profile-  see  malik  documents 
if  (nnp  .gt.  npt)  stop  ’too  many  points  in  profile’ 
write  (*,*)  ’working  profile  station  ’,nz,’  with  ’,nnp,’pts’ 
do  300  i  =  1 , iprnt  Inow  find  matching  prnt  station: 
if  (abs (sl-s (i) ) /si  .It.  small)  then 
jprnt  =  i 
go  to  301 
end  if 


230 

240 

250 
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300 


continue 

write  (*,*)  ’sl=  ’ ,sl 
stop  ’no  match  found’ 

301  continue 

write  (*,*)  ’found  match  at  print  station  ’,jprnt 
write  (*,*)  ’sl=  ’,sl,’  s(jprnt)=  ’,s(jprnt) 
resl  =  res(jprnt) 
rey  =  sqrt(resl) 
dstz  =  sl/rey 

* 

*  following  uses  radii  of  curvature 

* 

*  following  block  interpolates  contour  derivatives  for  radius  of  curvature 

*  from  original  sivells  data  file,  using  derivs  output  by  sivells 

if  (z(jprnt)  .gt.  0.0  .and.  iwherec  .eq.  1)  then 

*  aredownstream  of  throat,  so  do  radcurvature 

*  note  that  concave  curvature  should  have  a  minus  sign  for  Malik!! 

*  following  Numerical  Recipes  routine  finds  position  of  pt  in  array 

*  ssiv  that  is  just  below  si,  returns  in  jsiv 

call  locate(ssiv,nsiv,sl , jsiv) 

*  following  routine  performs  norder-pt  polynomial  interpolation 

call  polint (ssiv( j siv-norder+1) ,radcur (jsiv-norder+1) ,norder, 

>  si ,radcurl , errest) 

write  (*,*)  ’radcurtab=  ’ ,radcur(jsiv-norder+l) , 

>  ’  radcurinterpolated=  ’.radcurl 

if  (radcurl  .It.  huge  .and.  abs (errest) /radcurl  .gt.  eps)  then 
write  (*,*)  ’radcur  interpol.  err  est=  ’, errest 
write  (*,*)  ’  when  radcur=  ’, radcurl 
pause  ’  too  large?  ’ 
end  if 

epsxr  =  dstz/radcurl 

call  polint (ssiv(j siv-norder+1) , dydx(j siv-norder+1) ,norder, 

>  si, dydxl, errest) 

if  (abs (errest) /dydxl  .gt.  eps)  then 

write  (*,*)  ’dydx  interpol.  err  est=  ’, errest 
write  (*,*)  ’  when  dydx=  ’, dydxl 
pause  ’  too  large?  ’ 
end  if 

drdx  =  dydxl 
rmil  =  rmi(jprnt) 

else  if  (iwherec  .eq.  1)  then  ! upstream  of  throat  in  hopkins-hill  region 
epsxr  =  0.0  ! neglect  gortler  upstream  of  throat 

drdx  =  (rmi(jprnt+l)-rmi(jprnt))/(z(jprnt+l)-z(jprnt)) 
rmil  =  rmi(jprnt) 
radcurl  =  0.0  !flag 

else  ! iwherec  .eq.  2,  not  sivells,  use  original  data 
radcurl  =  radcur (jprnt)  ! local  value 

epsxr  =  dstz/radcurl 
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drdx  =  dydx(jprnt) 

rmil  =  rrmi(jprnt)  !use  value  interpolated  from  spline  fit 
end  if 

write  (3,350)  nz,z(jprnt) ,sl,rmi(jprnt) ,radcurl,epsxr 
350  f ormat (lx, i4,5(lx, lp,el3 . 6) ) 
c  rmil  =  rmi(jprnt)  !use  spline  fit  or  not,  depends,  moved  up 

thetal  =  theta(jprnt) 
del995  =  ye(jprnt) 
retlieta  =  res(jprnt)*thetal/sl 
amel  =  ame(jprnt) 
if  (amel  .gt.  0.)  then 
rethm  =  retheta/amel 
else 

write  (*,*)  ’amel=  ’  , amel ,  ’  nz=  1  ,nz 
stop  ’fatal  error’ 
end  if 

tel  =  te(jprnt) 
amuel  =  amue(jprnt) 
uel  =  ue(jprnt) 
xc  =  si 

pel  =  pe(jprnt) 

kodunit  =  0  (british  units 

igas  =  0  (perfect 

displc  =  dltast ( jprnt)  (displacement  thickness 

*  now,  write  general  info  to  file 

write  (2)  nz ,nnp,dstz ,rey, res 1 , epsxr, drdx, rmil , thetal , del995 , 

>  retheta, rethm, prandtl , kodunit , igas 
write  (2)  tel , amel , uel ,xc 

*  now  test  gortler  number  computations  vs.  emalik  style 

if  (epsxr  .It.  0.0)  then 

gortl  =  rey*sqrt (abs (epsxr) ) 
gortth  =  gortl* (thetal/dstz) **1 . 5 
else 

gortl  =  0.0 
gortth  =  0.0 
end  if 

write  (7,219)  12 . 0*z (jprnt) , si ,radcurl , gortl .gortth 

*  now,  read  the  profile  info: 

*  and  at  the  same  time  normalize 

xscal  =  del995*sqrt (resl) /si  (see  maliks  version  of  the  harris  code 
escall  =  (resl*amuel) / (rref *uref *sl*sqrt (2 . 0*xi ( jprnt) ) )  * 

>  rmil**jgeom  (rref  and  uref  added  8-30-91 
phi  =  atan(drdx)  (changed  9-8-94  sps 

*  yescal  changes  d/dytilde  derivs  to  d/dy/ye  derivs,  see  (15) 

yescal  =  del995/omega 

if  (nz  .eq.  ipro  .or.  nz  .eq.  1)  then 

write  (4,*)  ’debug  data  for  station  nz=  ’ ,nz 

write  (4,*)  ’resl , amuel , si ,xi=’ , resl , amuel , si ,xi (jprnt) 
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write  (4,*)  ’del995, xscal, escall=  ’ ,del995, xscal, escall 
write  (4,*)  ’phi,yescal  =  ’,phi,yescal 
write  (4,*)  ’rref,uref=  ’,rref,uref 
end  if 

if  (rmil  .le.  0.0)  then 

write  (*,*)  ’bltostab  debug:  rmil=  ’ ,rmil 

write  (*,*)  1  at  nz=  ’ ,nz, ’  sl=  ’.si,5  z(jprnt)=  ’,z(jprnt) 
write  (*,*)  1  rmi(jprnt)=  ’ ,rmi(jprnt) , 

>  ’  rrmi(jprnt)=  1 ,rrmi(jprnt) 

end  if 

do  400  j  =  l,nnp 

read  (1,*)  x(j) ,u(j) ,t(j) ,ul(j) ,tl(j)  Ireally  FZ  and  TZ,  accept  error 

*  following  rescales  eta  derivs  to  y/ye  derivs  for  malik  code 

*  this  next  pair  is  from  23b  in  harris  manual 

t23b(j)  =  1.0  +  x(j)*ye(jprnt)*cos(phi)/rmil 

*  this  next  line  derived  from  eqn  24b  in  manual,  changes  eta  derivs  to 

*  y/ye  derivs 

etascal(j)  =  escall  *  t23b(j)**jgeom  /  t(j) 

*  the  xscal  factors  in  the  following  are  to  convert  to  malik  code  form 

x(j)  =  x(j)*xscal  ! scaling  for  malik  code 
ul(j)  =  yescal+etascal (j ) *ul (j ) /xscal 
tl(j)  =  yescal+etascal (j ) *tl (j ) /xscal 
400  continue 

write  (*,*)  ’getting  derivatives’ 

call  scond(x,ul,u2,nnp)  !get  second  derivs  from  first 
call  scond(x,tl ,t2,nnp) 

*  ^change  to  use  of  utder  as  malik,  adapted  from  maliks 

*  call  utder (nnp,x,u,t ,ul ,u2 ,tl ,t2) 

*  don’t  have  derivs  from  harris  for  first  point 

*  now,  write  the  profile  info 

write  (2)  (x( j ) , j=l ,nnp) 
write  (2)  (u( j ) , j=l ,nnp) 
write  (2)  (ul(j) , j=l,nnp) 
write  (2)  (u2(j) , j=l,nnp) 
write  (2)  (t ( j ) , j=l ,nnp) 
write  (2)  (tl(j) , j=l,nnp) 
write  (2)  (t 2 ( j ) , j=l,nnp) 

*  for  checking 

ultest  =  (u(2)-u(l))/(x(2)-x(l)) 
if  (abs((ultest-ul(l))/ultest)  .gt.  eps)  then 

write  (*,*)  ’nz=  ’ ,nz, ’  ultest ,ul (1)=  ’, ultest ,ul (1) 
write  (*,*)  ’problems  with  generation  of  u  derivatives’ 
pause  ’looks  like  fatal  error’ 
end  if 

*  change  test  specs  due  to  profiles  being  so  flat,  adiabatic  wall  effects 

nnptest  =  nnp/2.0 

tltest  =  (t(nnptest)-t(nnptest-l))/(x(nnptest)-x(nnptest-l)) 
if  (tltest  .ne.  0.0)  then 
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if  ( (abs ( (tltest-tl (nnptest) ) /tltest)  .gt.  30*eps) 

>  .and.  (abs (tltest-tl (nnptest) )  .gt.  small))  then 


write 

(*,*) 

’nz=  ’,nz,’  tltest ,tl (nnptest)= 

> 

tltest ,tl (nnptest) 

write 

(*,*) 

’nnptest  =  ’ .nnptest 

write 

(*,*) 

’problems  with  generation  of  t  derivatives 

pause 

J looks 

like  fatal  error’ 

end  if 

end  if 

for  debug: 

if  (nz  .eq 

.  ipro 

.or.  nz  .eq.  1)  then 

call  scond(x,u,uldeb,nnp) 
write  (4,*)  ’nz=  1 ,nz, ’  ,  debug  info’ 
write  (4,*)  ’del995  (ye)=  ’,del995 
write  (4,*) 

>  ’x,u,ul ,uldeb ,uratio , t23b , etascal ,u2 ,t , tl , t2=  ’ , 

>  1 (as  written  to  bfl  file) 1 
do  900  i  =  l,nnp 

if  (uldeb(i)  .ne.  0.0)  then 
uratio  =  ul(i)/uldeb(i) 
else 

uratio  =  0.0  (arbitrary 
end  if 

write(4,850)  x(i) ,u(i) ,ul(i) ,uldeb(i) , uratio ,t23b(i) , 

>  etascal(i) ,u2(i) ,t(i) ,tl(i) ,t2(i) 

850  format (ll(lx,el8.12)) 

900  continue 
end  if 
1000  continue 

close (unit=4) 
close (unit=3) 

* 

stop 

end 

* - 

*  this  is  a  program  taken  from  sivells  to  compute  derivatives 

*  modified  7-13-90  to  deal  with  errors  in  endpoint 

SUBROUTINE  SC0ND  (A,B,C,KING) 

C  TO  OBTAIN  PARABOLIC  DERIVATIVE  OF  CURVE  (UNEQUALLY  SPACED  POINTS) 

*  IMPLICIT  REAL*8 (A-H , 0-Z) 

DIMENSION  A(*),  B(*),  C(*) 
data  eps/0.01/ 

N=KING-1 
DO  1  K=2 , N 

c  write  (*,*)  1 a( 1 ,k, 1 )=’ , a(k) 

S=A(K)-A(K-1) 

T=A(K+1)-A(K) 

1  C(K)=((B (K+l) -B (K))*S*S+(B(K)-B(K-1))*T*T)/(S*S*T+S*T*T) 
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SO=A (2) —A ( 1 ) 

if  (so  .eq.  0.)  then 

write  (*,*)  ’a(l ,2=’ ,a(l) ,a(2) 
stop  ’  SCOND :  so=0’ 
end  if 

T0=A(3)-A(2) 

if  (to  .eq.  0.)  stop  ’  SCOND:  to=0’ 

Q0=S0+T0 

C(1)=(-T0*(Q0+S0)*B(1)+Q0*Q0*B(2)-S0*S0*B(3) )/Q0/S0/T0 
*  following  added  when  got  bad  values  near  wall  sps  7-90 
clinear  =  (b(2)-b(l))/so 
if  (clinear  .ne.  0.0)  then 

error  =  abs(c(l)-clinear)/clinear 
else 

if  (c(l)  .ne.  0.0)  error  =  1.0 
end  if 

if  (error  .gt.  eps)  then 

write  (*,*)  ’SCOND:  problems  with  c(l)’ 
write  (*,*)  ’c(l) ,clinear=  ’, c (1) , clinear 
write  (*,*)  ’using  clinear’ 
c(l)  =  clinear 
end  if 

SF=A (KING-1) -A (KING-2) 

TF=A (KING) -A (KING-1 ) 

QF=SF+TF 

QST=QF*SF*TF 

C (KING) =  (SF* (QF+TF) *B (KING) -QF*QF*B (KING-1) +TF*TF*B (KING-2) ) /QST 

RETURN 

END 


* - 

*  The  subroutines  P0LINT  and  LOCATE  were  obtained  from  Numerical 

*  Recipes  by  Press  et .  al . ,  1st  edition. 

*  - 

*  following  routine  is  used  by  Frank  Chen’s  differentiation 

*  code  for  getting  curvatures  from  nozzle  contour  for  Gortler  work. 

*  this  code  put  into  the  BLT0STAB.F0R  program  sps  10-24-94 

* 

*  From  jerrypla@eagle.larc.nasa.gov  Thu  Oct  20  13:37  EST  1994 

*  Subject:  NOS  CSDS  CODE 

SUBROUTINE  CSDS (MAX , IX , X , F , DF , S , IPT , COEF , WK , I ERR) 


C*  * 


C*  PURPOSE:  * 
C*  SUBROUTINE  CSDS  FITS  A  SMOOTH  CUBIC  SPLINE  TO  A  * 
C*  UNIVARIATE  FUNCTION.  DATA  MAY  BE  UNEQUALLY  SPACED.  * 
C*  * 


C  E3.1 


C* 


* 
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c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


USE:  * 

CALL  CSDS (MAX , IX , X , F , DF , S , IPT , CDEF , WK , IERR)  * 

* 

MAX  INPUT  INTEGER  SPECIFYING  THE  MAXIMUM  NUMBER  OF  DATA  * 
POINTS  FOR  THE  INDEPENDENT  VARIABLE.  * 

* 

IX  INPUT  INTEGER  SPECIFYING  THE  ACTUAL  NUMBER  OF  DATA  * 
POINTS  FOR  THE  INDEPENDENT  VARIABLE.  IX#MAX.  * 

* 

X  ONE-DIMENSIONAL  INPUT  ARRAY  DIMENSIONED  AT  LEAST  * 

IX  IN  THE  CALLING  PROGRAM.  UPON  ENTRY  TO  CSDS,  * 

X(I)  MUST  CONTAIN  THE  VALUE  OF  THE  INDEPENDENT  * 

VARIABLE  AT  POINT  I.  * 

* 

F  ONE-DIMENSIONAL  INPUT  ARRAY  DIMENSIONED  AT  LEAST  * 

IX  IN  THE  CALLING  PROGRAM.  UPON  ENTRY  TO  CSDS,  * 

F (I)  MUST  CONTAIN  THE  VALUE  OF  THE  FUNCTION  AT  * 

POINT  X(I) .  * 

* 

DF  ONE-DIMENSIONAL  INPUT  ARRAY  DIMENSIONED  AT  LEAST  * 

IX  IN  THE  CALLING  PROGRAM.  UPON  ENTRY  TO  CSDS,  * 
DF(I)  MUST  CONTAIN  AN  ESTIMATE  OF  THE  STANDARD  * 

DEVIATION  OF  F(I) .  * 

* 

S  A  NON-NEGATIVE  INPUT  PARAMETER  WHICH  CONTROLS  THE  * 

EXTENT  OF  SMOOTHING.  S  SHOULD  BE  IN  THE  RANGE  * 

(IX- (2*IX) ** . 5) #S# (IX+ (2*IX) ** . 5) .  * 

* 

IPT  INPUT  INITIALIZATION  PARAMETER.  THE  USER  MUST  * 

SPECIFY  IPT=-1  WHENEVER  A  NEW  X  ARRAY  IS  * 

INPUT.  THE  ROUTINE  WILL  THEN  CHECK  TO  INSURE  THAT  * 
THE  X  ARRAY  IS  IN  STRICTLY  INCREASING  ORDER.  * 

* 

COEF  A  TWO-DIMENSIONAL  OUTPUT  ARRAY  DIMENSIONED  (MAX, 4)  * 

IN  THE  CALLING  PROGRAM.  UPON  RETURN,  COEF (I, J)  * 

CONTAINS  THE  J-TH  COEFFICIENT  OF  THE  SPLINE  FOR  * 

THE  INTERVAL  BEGINNING  AT  POINT  X(I) .  THE  * 

FUNCTIONAL  VALUE  OF  THE  SPLINE  AT  ABSCISSA  XI,  * 

WHERE  X(I)  . LE.  XI  .LE.  X(I+1) ,  IS  GIVEN  BY:  * 

F (XI ) = ( (COEF (1,4) *H+COEF (1,3)) *H+COEF (1,2)) *H  * 

+COEF(I , 1)  * 

WHERE  H=X1-X(I)  * 

* 

WK  A  ONE-DIMENSIONAL  WORK  AREA  ARRAY  DIMENSIONED  AT  * 

LEAST  (7*IX+9)  IN  THE  CALLING  PROGRAM.  * 

* 

IERR  OUTPUT  ERROR  PARAMETER:  * 

=0  NORMAL  RETURN.  NO  ERROR  DETECTED.  * 
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c* 

=J 

THE  J-TH  ELEMENT  OF  THE  X  ARRAY  IS  NOT  IN 

* 

c* 

STRICTLY  INCREASING  ORDER. 

* 

c* 

=-l 

THERE  ARE  LESS  THAN  FOUR  VALUES  IN  THE  X  ARRAY, 

.  * 

c* 

* 

c* 

UPON 

RETURN  FROM  CSDS,  THIS  PARAMETER  SHOULD  BE 

* 

c* 

TESTED  IN  THE  CALLING  PROGRAM. 

* 

c* 

* 

c* 

* 

c* 

* 

c* 

REQUIRED  ROUTINES 

-NONE 

* 

c* 

* 

c* 

LANGUAGE 

-FORTRAN 

* 

c* 

* 

c* 

DATE  RELEASED 

SEPTEMBER  5,  1973 

* 

c* 

* 

c* 

LATEST  REVISION 

MARCH  1975 

* 

c* ********************************************************************** 

c 

c 

DIMENSION  X(*) ,F(*) ,DF(*) ,COEF(MAX,*) ,WK(*) 

C 

C  SET  UP  WORKING  AREAS 

C 

IERR=0 

IF  (IPT  .NE.  -1)  GO  TO  5 

1  IPT=0 

IF(  IX  .LT.  4  )  GO  TO  2 
GO  TO  3 

2  IERR=-1 
RETURN 

3  1X1  =  IX-1 

DO  4  I  =  1,1X1 

IF  (  X(I  +  1)  -X(I)  .GT.  0  )  GO  TO  4 

I ERR  =  1+1 

RETURN 

4  CONTINUE 
NP1=IX  +1 
IB1  =  NP1 
IB2  =  IB1+NP1 
IB3  =  IB2+NP1+1 
IB4  =  IB3+NP1 
IB5  =  IB4+NP1 
IB6  =  IB5+NP1+1 
WK(1)  =  0. 

WK(2)  =  0. 

WK(IB2)  =  0. 

WK(IB3)  =  0. 

IJK2  =  IB2+NP1 
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WK(I JK2)=0 . 

I JK5  =  IB5  +  1 
WK(I JK5)=0 . 

I JK5  =  IB5  +  2 
WK(I JK5)=0 . 

WK(IB6)  =  0. 

I JK5  =  IB5+NP1 
WK(I JK5)=0 . 

5  CONTINUE 
P=0. 

H=X(2) -X (1) 

F2  =  -S 

FF=(F(2)-F(1))/H 
IF  (IX.LT.3)  GO  TO  25 
DO  6  1=3, IX 
G=H 

H=X (I) -X(I-l) 

E=FF 

FF=(F(I)-F(I-1))/H 
C0EF(I-1, 1)=FF-E 
I JK3  =  IB3+I 

WK (I JK3) = (G+H) * . 66666666666667 
I JK4  =  IB4+I 
WK (I JK4) =H/3 . 

I JK2  =  IB2+I 
WK (I JK2) =DF (1-2) /G 
WK(I)=DF(I)/H 
IJK1  =  IB1+I 

WK(IJK1)=-DF(I-1)/G-DF(I-1)/H 

6  CONTINUE 

DO  7  1=3, IX 
I JK1=IB1+I 
I JK2=IB2+I 

C0EF(I-1 , 2) =WK (I) *WK(I) +WK (I JK1) *WK(I JK1)+WK(I JK2) *WK (I JK2) 
C0EF(I-1 , 3) =WK (I) *WK(I JK1+1) +WK(I JK1) *WK(I JK2+1) 

C0EF(I-1 , 4) =WK (I) *WK(I JK2+2) 

7  CONTINUE 
C 

C  NEXT  ITERATION 

C 

10  IF  (IX.LT.3)  GO  TO  25 
DO  15  1=3, IX 

I JK1  =  IB1+I-1 
IJKO  =  1-1 

WK(IJK1)=FF*  WK(IJKO) 

I JK2  =  IB2+I-2 

IJKO  =  1-2 

WK (I JK2) =G*WK (I JKO) 
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I JKO  =  I 
I JK3  =  IB3+I 

WK (I JKO) =1 . / (P*COEF (1-1 , 2) +WK (I JK3) -FF*WK (I JK1) -G*WK(I JK2) ) 

I JK5  =  IB5+I 
I JKN  =  IJK5-1 
I JKO  =  IJKN-1 

WK (I JK5) =COEF (1-1 , 1) -WK(I JK1)  *  WK (I JKN) -WK (I JK2)  *WK(IJKO) 

I JK4  =  IB4+I 

FF=P*  COEF (1-1 , 3)+WK(I JK4) -H*  WK(IJKl) 

G=H 

H=COEF (1-1 , 4) *  P 
15  CONTINUE 
DO  20  1=3, IX 
J=IX-I+3 
IJK5  =  IB5+J 
I JK6  =  IJK5+1 
I JK7  =  IJK6+1 
IJK1  =  IB1+J 
I JK2  =  IB2+J 

WK(IJK5)  =  WK(J)*WK(IJK5)-WK(IJK1)*WK(IJK6)-WK(IJK2)*WK(IJK7) 
20  CONTINUE 
25  E=0 
H=0 
C 

C  COMPUTE  U  AND  ACCUMULATE  E 

C 

DO  30  1=2, IX 
G=H 

IJK5  =  IB5+I 

H  =  (WK(I JK5+1) -WK (I JK5) )/ (X(I)-X(I-l)) 

I JK6  =  IB6+I 

WK (I JK6) = (H-G) *  DF(I-l)  *  DF(I-l) 

E=E+WK (I JK6) * (H-G) 

30  CONTINUE 

G=-H*  DF (IX) *  DF (IX) 

IJK6  =  IB6+NP1 
WK(I JK6)=G 
E  =  E-G*H 
G=F2 

F2=E*P*P 

IF (F2 . GE . S  .OR.  F2.LE.G)  GO  TO  45 
FF=0 . 

IJK6  =  IB6+2 

H  =  (WK(I JK6+1) -WK(I JK6) )/ (X(2)-X(l)) 

IF  (IX  .LT.  3)  GO  TO  40 
DO  35  1=3, IX 
G=H 

I JK6  =  IB6+I 


54 


H  =  (WK(I JK6+1) -WK (I JK6) )/ (X(I)-X(I-l)) 

I JK1  =  IB1+I-1 
I JK2  =  IB2+I-2 

G  =  H-G-WK (I JK1) *WK (1-1) -WK (I JK2) *WK(I-2) 
FF=FF  +G  *  WK(I) *G 


WK (I)  = 
35  CONTINUE 
40  H=E-P*FF 
IF(H.LE.O) 


GO  TO  45 


C 

C 

C 

C 


C 

C 

C 

C 


UPDATE  THE  LAGRANGE  MULTIPLIER  P 
FOR  THE  NEXT  ITERATION 


P=P+ (S-F2) / ( (SQRT (S/E) +P) *H) 
GO  TO  10 


IF  E  LESS  THAN  OR  EQUAL  TO  S, 

COMPUTE  THE  COEFFICIENTS  AND  RETURN. 


45  DO  50  1=2, NP1 

I JK6  =  IB6+I 

COEFQ-1 , 1)  =F  (1-1)  -P*WK(I  JK6) 

I JK5  =  IB5+I 
COEF (1-1 , 3) =WK(I JK5) 

50  CONTINUE 

DO  55  1=2, IX 

H=X(I) -X(I-l) 

C0EF(I-1 , 4) = (COEF (I , 3) -C0EF(I-1 ,3) ) / (3 .  *H) 

C0EF(I-1 , 2) = (COEF (I , 1) -COEF (I- 1 , 1) ) /H  - (H*COEF (1-1 , 4)  +  COEF 
1  (1-1,3))  *  H 

55  CONTINUE 
9005  RETURN 
END 


55 


